\section{Introduction}
\label{sec:intro}

The Minkowski functional (MF) is one of the fundamental concepts in integral and stochastic geometry.
This is a series of geometric quantities defined for a bounded set in Euclidean space.
In the 2-dimensional case, the Minkowski functional of the set $M$ is a triplet consisting of the area $\mathrm{Vol}_2(M)$, the length of the boundary $\mathrm{Vol}_1(\partial M)$, and the Euler characteristic $\chi(M)$.
The Minkowski functional measures morphological future of $M$ in a different way from conventional moment-type statistics, and has been used in various scientific fields.

In cosmology, the Minkowski functional was introduced around in the 1990's,
and was first used to analyze the large-scale structure of the universe, and then the cosmic microwave background (CMB).
In particular, the Minkowski functional for the excursion set of the smoothed CMB map was firstly analyzed by \cite{schmalzing-buchert:1997} (cf.\,\cite{schmalzing-gorski:1998}).
The CMB radiation provides rich information on the very early stages of the universe.
Its signal is recognized as an isotropic and nearly Gaussian random field.
However, hundreds of inflationary models are available to deduce the various types of non-Gaussianity.
The Minkowski functional is used for the selection of such candidate models.

More precisely, let $X(t)$, $t\in E\subset\mathbb{R}^n$, be such a random field.
Its sup-level set with a threshold $v$,
\[
 E_v = \{t\in E \,|\, X(t)\ge v\} = X^{-1}([v,\infty)),
\]
is referred to as the \textit{excursion set} (Figure \ref{fig:excursion}).
Then, we can calculate the Minkowski functional curves $\{\mathcal{M}_j(E_v)\}_{j=0,\ldots,n}$ as a function of $v$.
When $X(t)$ is Gaussian, the expected Minkowski functional density (i.e., the MF per unit volume) has been well known to be described in the Hermite polynomial  (e.g., \cite{tomita:1986}).
Based on this fact, the departure between the sample Minkowski functional and the expected Minkowski functional is used as a measure of non-Gaussianity.

\begin{figure}
\begin{tabular}{ccc}
\scalebox{0.4}{\includegraphics{./rf-reduced.pdf}}
&\quad
\scalebox{0.4}{\includegraphics{./binary-reduced.pdf}}
&\quad
\scalebox{0.4}{\includegraphics{./triangle.pdf}}
\end{tabular}
\caption{A random field (left), excursion set (middle), and triangularization (right).}
\label{fig:excursion}
\end{figure}

On the other hand, the expected Minkowski functional for a non-Gaussian random field is unknown (except for the Gaussian related fields).
The non-Gaussianity of a random field is measured by the cumulant generating functional
\begin{align*}
\log &\,\mathbb{E}\Bigl[\exp\Bigl(\int_E \theta(t) X(t) \dd t\Bigr)\Bigr] \\
 &= \sum_{N}\frac{1}{N!} \int_{E^N} \theta(t_1)\cdots\theta(t_N)
 \mathrm{cum}(X(t_1),\ldots,X(t_N)) \dd t_1\ldots \dd t_N.
\end{align*}
To express weakly non-Gaussian random fields, we assume a weak non-Gaussianity through the $N$-point correlation function (the $N$-th cumulant)
\[
 \mathrm{cum}(X(t_1),\ldots,X(t_N)) = O\bigl(\nu^{N-2}\bigr), \quad N\ge 2,
\]
or equivalent conditions on higher-order spectra,
where $\nu\ll 1$ is a non-Gaussianity parameter.
In the context of CMB analysis, $\nu$ is the variation coefficient of the temperature measurement.
Such a weak non-Gaussian random field is realized by the weak non-linear filtering (e.g., the Wiener-chaos or Volterra series) of a Gaussian random field
 (see Section \ref{sec:numerical}).
Under this assumption,
\cite{matsubara:2003,hikage-komatsu-matsubara:2006} provided the perturbation formula for the expected Minkowski functional up to $O(\nu)$ including the skewness when the dimension is $n\le 3$,
 while \cite{matsubara:2010} provided the formula up to $O(\nu^2)$ including the skewness and kurtosis when $n=2$.
Then, they compared the theoretical Minkowski functionals under several candidate inflation models with observed and simulated CMB data to evaluate the goodness-of-fit.

The primary motivation of the paper is to complete the perturbation formulas when the dimension $n$ is arbitrary.
This not only completes the existing research, but also reveals the local power of the Minkowski functional for testing Gaussianity against contiguous alternatives.
Moreover, we shall characterize undetectable non-Gaussianities by the Minkowski functional statistic.
As previously stated, although the Minkowski functional is popular in morphological research, its performance as a test statistic has not yet to be studied.

In reality, the index set (i.e., survey area) $E$ is a bounded domain, and the boundary corrections should always be incorporated.
The Gaussian kinematic formula (GKF) established by \cite{taylor-adler:2009,adler-taylor:2007,adler-taylor:2019} provides the framework of the expected Minkowski functional on a bounded domain.
The GKF states that when the random field $X(t)$ is a Gaussian random field, or a Gaussian related field, the expected Minkowski functional is written as a linear combination of the Lipschitz--Killing curvatures $\mathcal{L}_j(E)$ of $E$ multiplied by the Euler characteristic density $\Xi_j(v)$ depending on the marginal distribution of $X(t)$.
Because of its simplicity, the GKF has been applied to cosmology recently (\cite{fantaye-marinucci:2015}, \cite{pranav-etal:2019:betti}).
In this paper, we show that the expected Minkowski functional for the weakly non-Gaussian random field has the same form as the GKF, although the weakly non-Gaussian is not a Gaussian related field even in the extended sense by \cite{panigrahi-etal:2019}.

The GKF was established in the research of the expected Euler characteristic heuristic,
which is the methodology of approximation for the tail probability of the maximum of a smooth random field.
The volume-of-tube method is almost the same method using a different base function.
In statistics, these methods and their variations are used in multiple comparisons, confidence band constructions, and change-point problems.
See, e.g.,
\cite{lu-kuriki:2017}, \cite{sommerfeld-sain-schwartzman:2018}, \cite{cheng-cammarota-cammarota-marinucci-schwartzman:2020} and references therein for recent developments.
In this paper, we will not pursue in this direction, but cite some papers from these areas for geometric preliminary.

The contributions of the paper can be summarized as follows:
(i) The expected Euler characteristic density for an isotropic field is provided as an asymptotic expansion up to $O(\nu^2)$ (Theorem \ref{thm:perturbation}).
(ii) The undetectable non-Gaussianity by the Minkowski functional is characterized (Theorem \ref{thm:loop}).
(iii) The perturbation formula for the expected Minkowski functional on a bounded domain is proved to have the same form as the GKF up to an arbitrary order (Theorem \ref{thm:boundary}).
(iv) As a crucial tool in perturbation analysis,
we pose identities on the Hermite polynomial (Lemmas \ref{lem:DThetadet1}, \ref{lem:DBdet}, and Corollary \ref{cor:Edet}).

The authors are preparing another cosmic paper \cite{matsubara-kuriki:2020},
where the perturbation formula for the Euler characteristic density  up to $O(\nu^2)$ is derived by an approach different from the one used in this paper.
The final formulas are presented in terms of the higher-order spectra for the purpose of convenience in cosmic research.

This paper is organized as follows.
Section \ref{sec:preliminary} is the preliminary part.
The expected Euler characteristic formula via Morse's theorem, and its simplified version under the isotropic assumption are briefly summarized (Sections \ref{sec:EEC} and \ref{sec:EEC-isotropic}).
Geometric notions and tools including the Minkowski functionals are introduced (Section \ref{sec:tube}).
Section \ref{sec:perturbation} covers the main results.
We first introduce the isotropic $N$-point correlation functions to represent the weak non-Gaussianity, and then discuss their properties (Sections \ref{sec:weak-gauss}, \ref{sec:derivatives}).
The perturbation formula for the Euler characteristic density up to $O(\nu^2)$ is derived (Section \ref{sec:EC_density}).
We also identify the undetectable derivatives of the $N$-point correlation functions for arbitrary $N$.
The expected Minkowski functional formula is proved to have the same form as the GKF (Section \ref{sec:boundary}).
We consider a weak non-Gaussian random field and determine the performance of the perturbation formula by a simulation study (Section \ref{sec:numerical}).
In Section \ref{sec:cosmic},
we discuss how the perturbation formula is used in the study of cosmology.
In the Appendix, the identities of the Hermite polynomial are proved (Section \ref{sec:Hermite}), and the regularity conditions are summarized (Section \ref{sec:regularity}).

\section{Preliminaries from integral/stochastic geometry}
\label{sec:preliminary}

Throughout the paper, we assume that $X(t)$ is a random field on $E\subset\mathbb{R}^n$ with mean 0, variance 1, and smooth sample path $t\mapsto X(t)$ in the following sense:
$X(t)$, $\nabla X(t)=(X_i(t))_{1\le i\le n}$, and $\nabla^2 X(t)=(X_{ij}(t))_{1\le i,j\le n}$ exist and are continuous with respect to $t=(t^1,\ldots,t^n)$ a.s., where
\[
 X_i(t) = \frac{\partial}{\partial t^i}X(t), \qquad
 X_{ij}(t) = \frac{\partial^2}{\partial t^i\partial t^j}X(t).
\]
In addition, we assume that $X(\cdot)$ is isotropic.
That is, the arbitrary finite marginal distribution $\{X(t)\}_{t\in E'}$, where $E'\subset\mathbb{R}^n$ is a finite set, is invariant under the group of rigid motion of $t$.

\subsection{Expected Euler characteristic via Morse's theory}
\label{sec:EEC}

In this subsection, we summarize the expected Euler characteristic formula without assuming that $X(\cdot)$ is isotropic.
Major materials were taken from \cite{adler-taylor:2007}, \cite{adler-taylor:2019}, \cite{takemura-kuriki:2002}, and \cite{schneider-weil:2008}.

For the index set $E$, we make two assumptions:

\begin{assumption}
\label{as:stratified}
 (i) $E$ is a $C^2$ piecewise smooth manifold \cite{takemura-kuriki:2002}, or equivalently, a $C^2$ stratified manifold \cite{adler-taylor:2007}.
Specifically, $E$ is disjointly divided as
\begin{equation*}
 E = \bigsqcup_{d=0}^{n}\partial_d E,
\end{equation*}
where $\partial_d E$ is a finite union of $d$-dimensional open $C^2$ submanifold.
$\partial_n E=\mathrm{int} E$ is the interior of $E$, and
$\partial_{n-1} E\sqcup\cdots\sqcup\partial_{0}E = \partial E$ is the boundary of $E$.

(ii) $E$ has a positive reach (positive critical radius).
\end{assumption}

The definition of the reach (critical radius) is provided in Section \ref{sec:tube}.
The formula for the reach can be found in \cite[Theorem 4.18]{federer:1959}.
See also Section 2.3 of \cite{kuriki-takemura:2009} for its derivation.
Roughly speaking, positive-reach means the local convexity.

For each point $t\in \partial_d E$, $E$ has a tangent cone $C_t=C_t(E)\subset\mathbb{R}^n$.
From the assumption of the positive-reach, $C_t(E)$ is always convex.
(The converse is not necessarily true.)
The dual of $C_t(E)$ in $\mathbb{R}^n$,
\[
 N_t = N_t(E) = \big\{ x\in \mathbb{R}^n \mid \langle x,y\rangle_{\mathbb{R}^n} \le 0,\ \forall y\in C_t(E) \bigr\},
\]
is referred to as the normal cone of $E$ at $t$.

Let $f$ be a real-valued $C^2$ function on $\mathbb{R}^n$.
The gradient and Hesse matrix of $f$ are denoted by $\nabla f$ and $\nabla^2 f$, respectively.
The restriction of $f$ on $\partial_d E$ is denoted by $f|_{\partial_d E}$.
The gradient and Hesse matrix of $f|_{\partial_d E}$ with respect to a standard coordinate $s=(s^1,\ldots,s^d)$ of $\partial_d E$ at $t$ are denoted by $\nabla f|_{\partial_d E}$ and $\nabla^2 f|_{\partial_d E}$, respectively.

The point $t=t^*\in\partial_d E$ such that $\nabla f|_{\partial_d E}(t^*)=0$ is called the critical point of $f|_{\partial_d E}$.

\begin{definition}
\label{def:morse}
$f$ is called Morse on $E$ if, at any critical point $t^*\in\partial_d E$ of $f|_{\partial_d E}$,

(i) $\nabla^2 f|_{\partial_d E}(t^*)$ is of full rank, and

(ii) $\nabla f(t^*)\ne 0$ if $d<n$.
\end{definition}

The following proposition is a version of Morse's theorem for a piecewise smooth manifold of positive-reach.
$\mathrm{sgn}(\cdot)$ is the sign function, and $\mathds{1}_{\{\cdot\}}$ is the indicator function.

\begin{proposition}[{\cite[Corollary 9.3.3]{adler-taylor:2007}}]
\label{prop:morse}
Suppose that $f$ is Morse on $E$.
Then,
\begin{equation}
\label{morse}
 \chi(f^{-1}((-\infty,v])) = \sum_{d=0}^n \sum_{t^*}\mathrm{sgn}\det\bigl(\nabla^2 f|_{\partial_d E}(t^*)\bigr)\mathds{1}_{\{f(t^*)\le v\}}\mathds{1}_{\{-\nabla f(t^*)\in N_t(E)\}},
\end{equation}
where the summation $\sum_{t^*}$ runs over the set of critical points of $f|_{\partial_d E}$.
\end{proposition}

We will use $-X$ as a Morse function under regularity conditions.
Using the Dirac measure $\delta(\cdot)$ on $\partial_d E$ with respect to the volume element $\dd s=\prod_{i=1}^d \dd s_i=\dd \mathrm{vol}_d(t)$ of $\partial_d E$, (\ref{morse}) with $f:=-X$ is formally rewritten as
\begin{align*}
 \chi(E_v)
&= \sum_{d=0}^n \int_{\partial_d E}\mathrm{sgn}\det(-\nabla^2 X|_{\partial_d E})\mathds{1}_{\{X(t)\ge v\}}\mathds{1}_{\{\nabla X(t)\in N_t\}} \delta(-\nabla X|_{\partial_d E}) |J_t| \dd s \\
&= \sum_{d=0}^n \int_{\partial_d E} \det(-\nabla^2 X|_{\partial_d E})\mathds{1}_{\{X(t)\ge v\}}\mathds{1}_{\{\nabla X(t)\in N_t\}} \delta(\nabla X|_{\partial_d E}) \dd \mathrm{vol}_d(t),
\end{align*}
where
$J_t = \det(\partial(-\nabla X|_{\partial_d E})/\partial s) = \det(-\nabla^2 X|_{\partial_d E})$
is the Jacobian.

Let $Q_t$ be an $n\times (n-d)$ matrix such that $Q_t Q_t^\top$ is the orthogonal projection onto $T_t(\partial_d E)^\perp$, the orthogonal complement of the tangent space of $\partial_d E$.
Then, $N_t \subset T_t^\perp(\partial_d E)$, and hence $\nabla f \in N_t$ $\Leftrightarrow$ $Q_t^\top\nabla f \in Q_t^\top N_t$.
By taking the expectation and changing the expectation and integral, we have
\begin{align}
\label{EEC}
 \mathbb{E}[\chi(E_v)]
=& \sum_{d=0}^n \int_{\partial_d E} \mathbb{E}[-\det(\nabla^2 X|_{\partial_d E})\mathds{1}_{\{Q_t \nabla X(t)\in Q_t N_t\}} \delta(\nabla X|_{\partial_d E})] \dd \mathrm{vol}_d(t) \\
=& \sum_{d=0}^n \int_{\partial_d E} \biggl[ \int_{\mathbb{R}^{1+n+\frac{1}{2}d(d+1)}} \det(-H)\mathds{1}_{\{x\ge v\}}\mathds{1}_{\{V_2\in Q_t^\top N_t\}} \nonumber \\
& \hfill \qquad\qquad\qquad \times
 p_t(x,0,V_2,H) \dd x \dd V_2 \dd H \biggr] \dd \mathrm{vol}_d(t), \nonumber
\end{align}
where $p_t(x,0,V_2,H)$ is the probability density of
$
 (X(t),\nabla X|_{\partial_d E}(t), Q_t^\top \nabla X(t), \allowbreak \nabla^2 X|_{\partial_d E}(t)) \in \mathbb{R}^{1+d+(n-d)+\frac{1}{2}d(d+1)}
$
for a fixed $t$.

The regularity conditions that assure $-X$ is a Morse function, and validate the integration including the Dirac delta function are given in \cite[Theorem 12.1.1]{adler-taylor:2007}, and are summarized in the Appendix.

\subsection{Expected Euler characteristic for isotropic random fields}
\label{sec:EEC-isotropic}

Thus far, we have not assumed that $X(t)$ is isotropic.
Hereafter, we assume that $X(t)$ is isotropic, and simplify formula  (\ref{EEC}) under the assumption.
We sometimes drop the argument $t$, because the marginal distribution is irrelevant to $t$.

The isotropic property implies that the marginal moment is irrelevant to $t$.
Without loss of generality, we let
\begin{equation}
\label{mean0var1}
 \mathbb{E}[X(t)]=0, \qquad \mathbb{E}[X(t)^2]=1.
\end{equation}
The covariance function is a function of the distance between the two points:
\begin{equation}
\label{rho}
 \mathbb{E}(X(t_1)X(t_2)) = \mathrm{cum}(X(t_1),X(t_2)) = \rho\bigl(\tfrac12\Vert t_1-t_2\Vert^2\bigr), \quad \rho(0)=1.
\end{equation}
To describe the smoothness of the covariance function, we let it be a function of the squared distance.
Under the regularity condition for exchanging the derivative $\partial/\partial t_i$ and the expectation $\mathbb{E}[\cdot]$, we have from (\ref{mean0var1}) that $\mathbb{E}[X_i]=0$, $\mathbb{E}[X_{ij}]=0$, $\mathbb{E}[X_i X]=0$, and $\mathbb{E}[X_i X_j]+\mathbb{E}[X X_{ij}]=0$
(see Section \ref{sec:regularity}).
Therefore, $\nabla X(t)$ is uncorrelated with $X(t)$ for each fixed $t$.
Moreover, from (\ref{rho}), we have
\[
 \mathbb{E}[X_i X_j] = \frac{\partial^2}{\partial t_{1}^{i}\partial t_{2}^{j}}\rho\bigl(\tfrac12\Vert t_1-t_2\Vert^2\bigr) \big|_{t_1=t_2} = -\rho'(0) \delta_{ij}.
\]
We change the variable from $\nabla^2 X$ to
\begin{equation*}
 R = (R_{ij})_{1\le i,j\le n} = \nabla^2 X + \gamma X I_n, \quad \gamma=-\rho'(0),
\end{equation*}
so that $R$ is uncorrelated to $X$ and $\nabla X$,
which simplifies the whole manipulation.

For a fixed $t\in\partial_d E$, let $P_t$ be an $n\times d$ matrix such that $P_t P_t^\top$ is the orthogonal projection matrix onto $T_t(\partial_d E)$.
Let $Q_t$ be an $n\times (n-d)$ matrix defined before.
Then, $P_t P_t^\top+Q_t Q_t^\top=I_n$.
We write the tangent cone and the normal cone as
\[
 C_t(E) = \mathbb{R}^d \oplus P_t C_{t,2}, \qquad N_t(E) = Q_t N_{t,2},
\]
respectively, where $C_{t,2}$ and $N_{t,2}$ are proper convex cones in $\mathbb{R}^{n-d}$.

Around $t\in\partial_d E$, we consider a geodesic coordinate system
$t^i=\phi^i(s)$, $s=(s^a)_{1\le a\le d}$ s.t.\ %
$\langle \partial\phi/\partial s^a, \partial\phi/\partial s^b\rangle |_t =\delta_{ab}$, and
$\langle \partial\phi/\partial s^a, \partial^2\phi/\partial s^b\partial s^c\rangle |_t = 0$.
Let $P_t=(\partial\phi^i/\partial s^a)_{1\le i\le n,1\le a\le d}$.

\begin{lemma}
\label{lem:2ndff}
Suppose that $t\in\partial_d E$ is a critical point of $X(t)$ restricted to $\partial_d E$.
That is,
 $V_1=P_t^\top\nabla X(t)=0$.
Then,
\[
 \nabla^2 X|_{\partial_d E}(t) =
 P_t^\top \nabla^2 X(t) P_t
 - S_t(V_2), \quad
 V_2 = Q_t^\top \nabla X(t),
\]
and $S_t(V_2)$ is the second fundamental form of $\partial_d E$ at $t$ with respect to the direction $V_2$.
\end{lemma}

\begin{proof}
When $t$ is a critical point,
\[
 \frac{\partial X(\phi(s))}{\partial s^a}
 = \sum_{k=1}^n X_k\frac{\partial \phi^k(s)}{\partial s^a} = (P_t\nabla X(t))_a =0
\]
for $1\le a\le d$, and
\begin{align*}
 (\nabla|_{\partial_d E}^2X(t))_{ab}
 = \frac{\partial^2 X(\phi(s))}{\partial s^a\partial s^b}
 =& \sum_{k,l=1}^n X_{kl}(t)\frac{\partial \phi^k(s)}{\partial s^a}\frac{\partial \phi^l(s)}{\partial s^b}
 + \sum_{k=1}^n X_k(t)\frac{\partial^2 \phi^k(s)}{\partial s^a\partial s^b} \\
 =& (P_t^\top \nabla^2 X(t) P_t)_{ab}
  + \Bigl\langle Q_t V_2,\frac{\partial^2 \phi(s)}{\partial s^a\partial s^b} \Bigr\rangle.
\end{align*}
\end{proof}

\begin{lemma}
\label{lem:Vbar}
Let $V_1$ and $V_2$ be defined in Lemma \ref{lem:2ndff}.
$\Bar V_{2}=V_{2}/\Vert V_{2}\Vert$ is distributed uniformly on the unit sphere of $Q_t^\top T_t(\partial E_d)^\perp$, independently from $(X,V_{1},\Vert V_2\Vert,R_1)$, where $R_1 = P_t^\top R P_t$.
\end{lemma}
\begin{proof}
The assumption of isotropic property implies the invariance of distribution under $(X,V_{1},V_{2},R_1)\mapsto (X,V_{1},P_2 V_{2},R_1)$, where $P_2\in O(n-d)$.
The conditional distribution of $V_2$ when $(X,V_1,R_1)$ is given is orthogonally invariant, and the conditional distribution of $\Bar V_2$ is a uniform distribution on the unit sphere of $Q_t^\top T_t(\partial E_d)^\perp$ independently from $(X,V_1,R_1)$ and $\Vert V_2\Vert$.
\end{proof}

By choosing an appropriate orthonormal coordinate again, we can assume that
\begin{align*}
 S_t(V_{2})
=& S_t\bigl(\Bar V_{2}\bigr)\Vert V_{2}\Vert
= \mathrm{diag}\bigl(\sigma_{t,1}\bigl(\Bar V_{2}\bigr),\ldots,\sigma_{t,k}(\Bar V_{2})\bigr)\Vert V_{2}\Vert,
\end{align*}
where $\sigma_{ti}(\Bar V_{2})$, $i=1,\dots,k$, are principal curvatures of $\partial_d E$ at $t$ with respect to the normal direction $Q_t \Bar V_{2}\in N_t(E)$.
For a subset $I=\{i_1,\ldots,i_k\}\subset\{1,\ldots,d\}$ ($i_1<\cdots<i_k$) such that $|I|=k$, let $\Bar I=\{1,\ldots,d\}\setminus I$, and
\[
 R_{\Bar I} = ((R_1)_{ij})_{i,j\in \Bar I}, \quad
 S_{t,I}\bigl(\Bar V_{2}\bigr) = \mathrm{diag}\bigl(\sigma_{t,i_1}\bigl(\Bar V_{2}\bigr),\ldots,\sigma_{t,i_k}(\Bar V_{2})\bigr).
\]
The determinant in (\ref{EEC}) is expanded as
\[
 \det\bigl(-R_1 +\gamma X I_d + S_t(V_{2})\bigr)
= \sum_{e=0}^d \sum_{I: |I|=e} \det(-R_{\Bar I}+\gamma X I_{d-e}) \det S_{t,I}\bigl(\Bar V_{2}\bigr)\Vert V_{2}\Vert^e.
\]
We substitute this into (\ref{EEC}),
and find that $\Bar V_2$ is uniformly distributed on the unit sphere $\S(\mathbb{R}^{n-d})$ of $\mathbb{R}^{n-d}$ independently from $X, R_{\Bar I}$, $V_1$ and $\Vert V_2\Vert$.
Let $R_{11}=((R_1)_{ij})_{1\le i,j\le d-e}$.
Then,
\begin{align}
 \mathbb{E}[\chi(E_v)] = \sum_{d=0}^n \sum_{e=0}^d L_{n-d,e} \Xi_{d,e}(v),
\label{LXi}
\end{align}
where
\begin{align}
\label{Xide}
 \Xi_{d,e}(v) =
 \int_{x=v}^\infty \biggl[ & \int_{\mathrm{Sym}(d-e)} \det\bigl(-R_{11}+\gamma x I_{d-e} \bigr) \\
 & \times \int_{\mathbb{R}^{n-d}} \Vert V_2\Vert^e p(x,(0,V_{2}),R_{11}) \dd V_{2} \dd R_{11} \biggr] \dd x, \nonumber
\end{align}
and
\begin{align}
 L_{n-d,e} =
 \int_{\partial E_d} \mathbb{E}\bigl[ \mathds{1}(\xi\in N_{t,2}) \mathrm{tr}_e S_t(\xi)\bigr] \dd \mathrm{vol}_d(t), \quad \xi\sim\mathrm{Unif}\bigl(\S(\mathbb{R}^{n-d})\bigr)
\label{Lde}
\end{align}
with
\[
  \mathrm{tr}_e S_t(\xi) = \sum_{I: I\subset\{1,\ldots,d\},\,|I|=e}\det S_{t,I}(\xi).
\]

When the leading term from $\mathrm{int}(E)$ is dominant and the contribution of the boundaries $\partial_d E$ ($d<n$) is negligible, (\ref{LXi}) is approximated as
\begin{equation*}
 \mathbb{E}[\chi(E_v)] \approx L_{n,0} \Xi_{n}(v),
\end{equation*}
where
\begin{equation}
\label{Xin}
 \Xi_{n}(v) =
 \int_{x=v}^\infty \biggl[ \int_{\mathrm{Sym}(n)} \det\bigl(-R +\gamma x I_n\bigr) p(x,0,R) \dd R \biggr] \dd x,
\end{equation}
$p(x,0,R)$ is the probability density of $(X,V,R)$, and
\[
 L_{n,0} = \int_E \dd t = \mathrm{Vol}_n(E).
\]
We call $\Xi_{n}(v)$ in (\ref{Xin}) as the Euler characteristic density.
In Section \ref{sec:perturbation}, we obtain $\Xi_{n}=\Xi_{n,0}$ in (\ref{Xin}), $\Xi_{d,e}$ in (\ref{Xide}), and $\mathbb{E}[\chi(E_v)]$ in (\ref{LXi}), in turn.
Because of the assumption of isotropicity, we can let $P_t^\top=(I_d,0)$, $Q_t^\top=(0,I_{n-d})$ included in the integrals (\ref{Xide}) and  (\ref{Xin}).

\subsection{Tube, critical radius, and Minkowski functional}
\label{sec:tube}

The critical radius is a geometric quantity for a subset $M$ in Euclidean space $\mathbb{R}^n$.
This is defined via the tube about $M$ below.

Tube about $M$ with radius $\rho$ is defined by a set of points whose distance from $M$ is less than or equal to $M$:
\begin{equation}
\label{tube}
 \mathrm{Tube}(M,\rho) = \bigl\{ x\in\mathbb{R}^n \mid \mathrm{dist}(x,M) \le \rho \bigr\}, \qquad
 \mathrm{dist}(x,M) = \min_{y\in M}\Vert y-x \Vert.
\end{equation}

Suppose that the radius $\rho$ increased from 0 to infinity.
Unless $M$ is convex, there exists a finite value of $\rho=\rho_c(M)$ such that the tube $\mathrm{Tube}(M,\rho)$ has a self-overlap.
This threshold is referred to as the {\it critical radius} or the {\it reach} of $M$:
\[
 \rho_c(M) = \inf\bigl\{\rho\ge 0 \mid \mathrm{Tube}(M,\rho) \mbox{ has a self-overlap} \bigr\}.
\]
$M$ is said to be of positive-reach if $\rho_c(M)$ is strictly positive.
The classical Steiner formula states that, when $M$ is a $C^2$ piecewise smooth manifold of positive-reach, for all $\rho\in[0,\rho_c(M)]$, the volume of the tube (\ref{tube}) is expressed as a polynomial in $\rho$:
\begin{equation}
\label{steiner}
 \mathrm{Vol}_n(\mathrm{Tube}(M,\rho))
 = \sum_{j=0}^n \omega_{n-j}\rho^{n-j} \mathcal{L}_j(M)
 = \sum_{j=0}^n \rho^{j} \binom{n}{j} \mathcal{M}_j(M),
\end{equation}
where
\begin{equation*}
 \omega_d=\mathrm{Vol}_d(\mathbb{B}^d)=\frac{\pi^{d/2}}{\Gamma(d/2+1)}
\end{equation*}
is the volume of the unit ball in $\mathbb{R}^d$.
The Minkowski functional
$\mathcal{M}_j(M)$ of $M$, and the Lipschitz--Killing curvature $\mathcal{L}_j(M)$ of $M$ are defined as the coefficients of the polynomial.
$\mathcal{L}_j(M)$ is defined independently from the dimension $n$ of the ambient space.
There are variations of the definitions of the Minkowski functional.
The definition in (\ref{steiner}) is by
\cite[Section 14.2]{schneider-weil:2008}.

Owing to the Gauss-Bonnet theorem,
the Minkowski functional of the largest degree is proportional to the Euler characteristic (EC) of $M$:
\[
 \chi(M) = \mathcal{L}_0(M) = \mathcal{M}_n(M)/\omega_n.
\]

The next two propositions are used in Section \ref{sec:boundary}.
Let
\begin{equation*}
 \Omega_d = \mathrm{Vol}_{d-1}(\S^{d-1}) = \frac{2\pi^{d/2}}{\Gamma(d/2)}
\end{equation*}
be the volume of the unit sphere in $\mathbb{R}^d$.
Note that $\omega_d=\Omega_d/d$.

\begin{proposition}[{Crofton's intersection formula, \cite[Theorem 5.1.1]{schneider-weil:2008}}
\label{prop:crofton}
Let $A(n,k)$ be the set of $k$-dimensional affine subspaces in $\mathbb{R}^n$.
Let $\mu_{n,k}$ be the invariant measure on $A(n,k)$ with respect to the isotropic group normalized by
\[
 \mu_{n,k}(\{L\in A(n,k) \,|\, L\cap \mathbb{B}^n\ne\emptyset \}) = \omega_{n-k}.
\]
Then,
\begin{equation*}
 \int_{A(n,n-k)} \mathcal{L}_j(M\cap L) \dd \mu_{n,n-k}(L) = c_{n,j,k} \mathcal{L}_{n+j-k}(M),
\end{equation*}
where
\begin{equation*}
 c_{n,j,k} = \frac{\Gamma\bigl(\frac{k+1}{2}\bigr)\Gamma\bigl(\frac{k+j+1}{2}\bigr)}{\Gamma\bigl(\frac{j+1}{2}\bigr)\Gamma\bigl(\frac{n+1}{2}\bigr)}.
\end{equation*}
\end{proposition}

\begin{proposition}
\label{prop:LK_curvature}
Let $E=\bigsqcup_{d=0}^n\partial_d E\subset\mathbb{R}^n$ be the $C^2$ piecewise smooth manifold of positive-reach.
Then, the Lipschitz--Killing curvature of $E$ is
\begin{align*}
 \mathcal{L}_d(E) = \sum_{e=d}^n L_{n-e,e-d} \frac{\Omega_{n-e}}{\Omega_{n-d}},
\end{align*}
where $L_{n-e,e-d}$ is defined in (\ref{Lde}).
\end{proposition}

\begin{proof}
According to the volume formula of tube,
\begin{equation*}
\mathrm{Vol}_n(\mathrm{Tube}(E,\rho))
= \sum_{d=0}^n \int_{\partial E_d} \int_{r<\rho} \int_{v\in \S(N_{t,2})} \det\bigl( I_d + r S_t(v)\bigr) r^{n-d-1} \dd v \dd r \dd t,
\end{equation*}
where $\dd v$ is the volume element of $\S(N_{t,2})=\{v/\Vert v\Vert \,|\, v\in N_{t,2}\}$.
Note that $\dd v/\Omega_{n-d}$ is the uniform distribution on $\S(N_{t,2})$,
$\int_{v\in \S(N_{t,2})} \mathrm{tr}_e(S_t(v)) \dd v = L_{n-d,e}\Omega_{n-d}$,
and
\begin{align*}
\mathrm{Vol}_n(\mathrm{Tube}(E,\rho))
=& \sum_{d=0}^n \sum_{e=0}^d L_{n-d,e} \Omega_{n-d} \int_{r<\rho} r^{n-d+e-1} \dd r \\
=& \sum_{d=0}^n \sum_{e=0}^d L_{n-d,e} \frac{\Omega_{n-d}}{\Omega_{n-d+e}} \times \omega_{n-d+e} \rho^{n-d+e} \\
=& \sum_{d=0}^n \Biggl(\sum_{e=d}^n L_{n-e,e-d} \frac{\Omega_{n-e}}{\Omega_{n-d}}\Biggr) \omega_{n-d} \rho^{n-d}
= \sum_{d=0}^n \mathcal{L}_d(E) \omega_{n-d} \rho^{n-d}.
\end{align*}
\end{proof}

\section{Perturbation of the Minkowski functional}
\label{sec:perturbation}

Throughout the paper, we assume that $X(t)\in\mathbb{R}$, $t\in E\subset\mathbb{R}^n$, is a zero-mean, unit-variance random field with continuous derivatives $\nabla X(t)$ and $\nabla^2 X(t)$ with probability one.
We also assume that $X(t)$ is isotropic (except for Section \ref{sec:EEC}).

\subsection{Weakly non-Gaussian random fields}
\label{sec:weak-gauss}

To describe the weak non-Gaussianity, we introduce the $N$-point correlation functions
defined as a joint cumulant of the $N$ variable $(X(t_1),\ldots,X(t_N))$, where $t_i$'s are arbitrary points.
In cosmology, the weakly nonlinear evolution of cosmic fields is modeled such that the $N$-point correlation is of the order $O(\nu^{N-2})$, where $\nu\ll 1$ is the variation coefficient of the cosmic fields (i.e., CMB, large-scale structure, and weak lensing fields)  \cite{matsubara:2003}.
This asymptotics is the same as the setting of the central limit theorem (CLT) for independent and identically distributed samples when $1/\nu^2$ is the sample size.
Therefore, we will discuss the perturbation theory in the framework of the CLT (\cite{petrov:1975,bhattacharya-Rao:2010}).

Because of the isotropic property, we assume the $N$-point correlation function:
\begin{align*}
 \mathrm{cum}(X(t_1),\ldots,X(t_N)) = \nu^{N-2} \kappa^{(N)}\bigl((\tfrac12\Vert t_a-t_b\Vert^2)_{1\le a<b\le N}\bigr),
\end{align*}
where $\kappa^{(N)}$ is a smooth function.
Its argument is an $\binom{N}{2}$-dimensional vector with elements arranged by the lexicographic order of $(a,b)$.
The 2-point correlation function $\kappa^{(2)}(\cdot)$ is the covariance function $\rho(\cdot)$.
The 3- and 4-point correlation functions are
\begin{align*}
\mathrm{cum}(X(t_1),X(t_2),X(t_3)) =& \mathbb{E}(X(t_1)X(t_2)X(t_2)) \\
 =& \nu \kappa^{(3)}\bigl(\tfrac12\Vert t_1-t_2\Vert^2,\tfrac12\Vert t_1-t_3\Vert^2,\tfrac12\Vert t_2-t_3\Vert^2\bigr),
\end{align*}
and
\begin{align*}
\mathrm{cum}( & X(t_1),X(t_2),X(t_3),X(t_4)) \\
& = \mathbb{E}(X(t_1)X(t_2)X(t_3)X(t_4)) -\mathbb{E}(X(t_1)X(t_2))\mathbb{E}(X(t_3)X(t_4))[3] \\
& = \nu^2 \kappa^{(4)}\bigl(\tfrac12\Vert t_1-t_2\Vert^2,\tfrac12\Vert t_1-t_3\Vert^2,\ldots,\tfrac12\Vert t_3-t_4\Vert^2\bigr),
\end{align*}
where ``$[3]$'' means that it represents three symmetric terms.
Recall that $\kappa^{(N)}$ is invariant with respect to the permutation $S_N$ of $t_i$'s.
As a result, $\kappa^{(3)}(x_{12},x_{13},x_{23})$ is symmetric in its arguments, but $\kappa^{(N)}$ ($N\ge 4$) are not symmetric.

To evaluate the Euler characteristic densities $\Xi_n$ in (\ref{Xin}) or $\Xi_{d,e}$ in (\ref{Xide}), we need to identify the joint distribution of $(X,\nabla X,\nabla^2 X)$ and then $(X,V,R)$ with $V=\nabla X$, and $R=\nabla^2 X+\gamma X I_n$.

Using the parameters
$T = (\tau_i)\in\mathbb{R}^n$ and $\Theta = (\theta_{ij})\in\mathrm{Sym}(n)$ such that
\begin{equation}
\label{Theta}
 \theta_{ij} = \frac{1+\delta_{ij}}{2}\tau_{ij} \quad (i\le j),
\end{equation}
the characteristic function of $1+n+n(n+1)/2=\binom{n+2}{2}$ dimensional random variables $(X,\nabla X,\nabla^2 X)$ is $\psi(i s,i T,i \Theta)$, where
\begin{align}
\label{mgf_X''}
\psi(s,T,\Theta)
 =& \mathbb{E}\Bigl[ e^{s X+\sum_i \tau_i X_i+\sum_{i\le j} \tau_{ij}X_{ij}} \Bigr]
 = \mathbb{E}\Bigl[ e^{s X+\langle T,\nabla X\rangle+\mathrm{tr}(\Theta \nabla^2 X)} \Bigr]
\end{align}
is the moment generating function (mgf).
Then, the mgf of $(X,\nabla X,R)$ is
\begin{equation*}
 \mathbb{E}\Bigl[ e^{s X+\langle T,\nabla X\rangle+\mathrm{tr}(\Theta \nabla^2 R)}\Bigr]
 = \psi(\s,T,\Theta), \quad \s=s+\gamma\mathrm{tr}(\Theta).
\end{equation*}

Under smoothness conditions, the cumulant generating function $\log\psi(s,T,\Theta)$ has a Taylor series
\[
 \log\psi(s,T,\Theta) = \sum_{u,v,w\ge 0}\frac{\nu^{u+v+w-2}}{u!\,v!\,w!} K_{u,v,w}(s,T,\Theta),
\]
where, for $N=u+v+w$,
\begin{align}
\label{Kuvw}
& \nu^{N-2} K_{u,v,w} (s,T,\Theta) \\
& = \mathrm{cum}\bigl(\underbrace{s X,\ldots,s X}_u,\underbrace{\langle T,\nabla X\rangle,\ldots,\langle T,\nabla X\rangle}_v,\underbrace{\mathrm{tr}(\Theta\nabla^2 X),\ldots,\mathrm{tr}(\Theta\nabla^2 X)}_w\bigr) \nonumber \\
& = \nu^{N-2} s^u \biggl(\prod_{b=u+1}^{u+v} \langle T,\nabla_{t_b}\rangle \prod_{c=u+v+1}^{N} \mathrm{tr}\bigl(\Theta\nabla^2_{t_c}\bigr)\biggr) \kappa^{(N)}\bigl((\tfrac{1}{2}\Vert t_a-t_b\Vert^2)_{a<b}\bigr)\Big|_{t_1=\cdots=t_{N}} \nonumber
\end{align}
with $\nabla_{t_b}=(\partial/\partial t_{b}^{i})_{1\le i\le n}$, $\nabla^2_{t_c}=(\partial^2/\partial t_{c}^{i}\partial t_{c}^{j})_{1\le i,j\le n}$.

The form of $K_{u,v,w}(s,T,\Theta)$ is specified by the isotropic assumption.

\begin{lemma}
\label{lem:invariant_poly}
$K_{u,2v,w}(s,T,\Theta)$ is a linear combination of
\begin{equation}
\label{K}
 s^u \prod_{j\ge 0}(T^\top \Theta^j T)^{v_j} \prod_{k\ge 1}\mathrm{tr}(\Theta^k)^{w_k},
\end{equation}
where $v_j$ and $w_k$ are non-negative integers satisfying
\[
 \sum_{j\ge 0}v_j = v, \quad \sum_{j\ge 1} j v_j + \sum_{k\ge 1} k w_k = w.
\]
\end{lemma}

\begin{proof}
By considering the orthogonal transformation $t\mapsto P t$, $P\in O(n)$, we have the invariance in distribution
 $(X,\nabla X,\nabla^2 X)\mathop{=}^d (X,P\nabla X,P\nabla^2 X P^\top)$.
Because of the multilinearity of the cumulant, $K_{u,v,w}(s,T,\Theta)=s^u K_{u,v,w}(1,T,\Theta)$ and $K_{u,v,w}(1,P^\top T,P^\top \Theta P) = K_{u,v,w}(1,T,\Theta)$.
By letting $P=-I_n$, we see that $K_{u,v,w}(1,T,\Theta)=0$ when $b$ is odd.
When $b$ is even, $K_{u,v,w}(1,T,\Theta)$ is a quadratic function of $T$, and can be regarded as a function of $(t_i t_j)_{1\le i,j\le n}= T T^\top\in\mathrm{Sym}(n)$.

As an extension of the zonal polynomial, \cite{davis:1980} introduced an invariant polynomial of two symmetric matrices $A,B\in\mathrm{Sym}(n)$ which is invariant under the transformation $(A,B)\mapsto (P^\top A P,P^\top B P)$.
$s^u K_{u,v,w}(1,T,\Theta)$ is an invariant polynomial in $(T T^\top,\Theta)$, which is a linear combination of (\ref{K}) as shown in (4.8) of \cite{davis:1980}.
\end{proof}

In the limit $\nu\downarrow 0$, uncorrelated random variables $X$, $\nabla X$, and $\nabla^2 X+\gamma I_n=R$ are independently distributed as the Gaussian distributions $X\sim N(0,1)$, $\nabla X\sim N_n(0,\gamma I_n)$, and $R\in\mathrm{Sym}(n)$ is a zero-mean Gaussian random matrix with the covariance structure
\begin{align*}
\mathrm{cum}(R_{ij},R_{kl})
= \alpha \frac{1}{2} (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) + \beta \delta_{ij}\delta_{kl},
\end{align*}
where
\begin{equation*}
 \alpha = 2\rho''(0), \quad \beta = \rho''(0)-\rho'(0)^2.
\end{equation*}
The probability density function (pdf) of the limiting $R$ with respect to $\dd R=\prod_{i\le j} \dd R_{ij}$ is proportional to
\begin{align*}
\exp\Bigl\{-\frac{1}{2\alpha}\mathrm{tr}\bigl((R-\tfrac{1}{n}\mathrm{tr}(R) I_n)^2\bigr) - \frac{1}{2 n(\alpha+n\beta)}\mathrm{tr}(R)^2 \Bigr\}.
\end{align*}
This is a proper density if
\begin{equation}
\label{range}
 \alpha>0,\,\alpha+n\beta>0, \quad\mbox{or equivalently},\quad
 \rho''(0) > \frac{n}{n+2}\rho'(0)^2 > 0.
\end{equation}
This random matrix $R$ is refereed to as the Gaussian orthogonal invariant matrix (\cite{cheng-schwartzman:2018}).

\subsection{Derivatives of point correlation functions}
\label{sec:derivatives}

Here, we will evaluate $K_{u,v,w}(s,T,\Theta)$ in detail, specifically up to $N=u+v+w\le 4$.
Lemma \ref{lem:invariant_poly} states these are linear combinations of (\ref{K}) with coefficients the partial derivatives of the point correlation function evaluated at the origin
\begin{equation}
\label{kappaN}
\kappa^{(N)}_{(i_{12},\ldots,i_{N-1,N})}(0,\ldots,0),
\end{equation}
where
\begin{align*}
 \kappa^{(N)}_{(i_{12},\ldots,i_{N-1,N})}(x_{12},\ldots,x_{N-1,N}) =
 \biggl(\prod_{1\le a<b\le N}\Bigl(\frac{\partial}{\partial x_{ab}}\Bigr)^{i_{ab}}\biggr) \kappa^{(N)}(x_{12},\ldots,x_{N-1,N}).
\end{align*}
In the following, we evaluate the cumulants up to the 4th order using the 3- and 4-point correlation functions $\kappa^{(3)}$ and $\kappa^{(4)}$, respectively.
Table \ref{tab:derivatives} lists the derivatives (\ref{kappaN}) that appear in the expression $K_{u,v,w}$ for $N=u+v+w\le 4$, as well as their abbreviated notations.
Because $\kappa^{(3)}(x_{12},x_{13},x_{23})$ is symmetric in its arguments, there are three 1st derivatives with the same value
$\kappa^{(3)}_{(1,0,0)}(0,0,0)=\kappa^{(3)}_{(0,1,0)}(0,0,0)=\kappa^{(3)}_{(0,0,1)}(0,0,0)$.
The number of such derivatives is referred to as the \textit{multiplicity}.

For the derivative (\ref{kappaN}), the \textit{diagram} is defined as an undirected graph with vertex set $\{1,\ldots,N\}$ and edge set
\[
 \{ \underbrace{(1,2),\ldots,(1,2)}_{i_{12}\,\rm times},\underbrace{(1,3),\ldots,(1,3)}_{i_{13}\,\rm times},\ldots\ldots \}.
\]
This edge set is a multiset, which means that multiple edges connecting a pair of two vertices are allowed.
The vertices of degree 0 (vertices without edges) can be omitted from the diagram.

The diagram is characterized as follows.
\begin{lemma}
\label{lem:diagram}
The diagram of (\ref{kappaN}) is an undirected graph with a maximum of $N$ vertices.
Each connected component is a tree, or graph containing one loop with a length greater than or equal to $2$.
\end{lemma}

\begin{proof}
Let $x_{ab} = \Vert t_a-t_b\Vert^2/2 = \sum_{i=1}^n (t_{a}^{i}-t_{b}^{i})^2/2$.
Since $\partial/\partial t_{a}^{i}
=\sum_{a'\ne a}(t_{a}^{i}-t_{a'}^{i})\partial/\partial x_{aa'}$,
we have $\partial/\partial t_{a}^{i}|_{t_1=\cdots=t_N}=0$, and
\[
 \partial^2/(\partial t_{a}^{i}\partial t_{b}^{j})|_{t_1=\cdots=t_N}=
 \delta_{ij}\times\begin{cases}
 \sum_{a'\ne a}\partial/\partial x_{aa'} & (a=b), \\
 - \partial/\partial x_{ab} & (a\ne b).
 \end{cases}
\]
From this observation,
we find that the derivative $\langle T,\nabla_{t_{b_1}}\rangle\langle T,\nabla_{t_{b_2}}\rangle$ yields the edge $(b_1,b_2)$, while the derivative $\mathrm{tr}(\Theta\nabla^2_{t_c})$ yields the edge $(c,c')$, where $c'$ is any vertex other than $c$.
The diagram is constructed by using these edges as the building blocks.
We see that a loop-free diagram consists of one $(b_1,b_2)$-type edges and $(c,c')$-type edges,
and a diagram containing a loop consists only of $(c,c')$-type edges.
(For example,
$\langle T,\nabla_{t_1}\rangle\langle T,\nabla_{t_4}\rangle\mathrm{tr}(\Theta\nabla^2_{t_2}) \allowbreak \mathrm{tr}(\Theta\nabla^2_{t_3})\mathrm{tr}(\Theta\nabla^2_{t_5})\mathrm{tr}(\Theta\nabla^2_{t_6})$ produces Figure \ref{fig:loop+tree} (left), and
$\mathrm{tr}(\Theta\nabla^2_{t_1})\mathrm{tr}(\Theta\nabla^2_{t_2}) \allowbreak \mathrm{tr}(\Theta\nabla^2_{t_3})\mathrm{tr}(\Theta\nabla^2_{t_4})\mathrm{tr}(\Theta\nabla^2_{t_5})\mathrm{tr}(\Theta\nabla^2_{t_6})$ produces Figure \ref{fig:loop+tree} (right).)
Hence, the resulting diagram is a union of connected graphs, and each connected graph has at most one loop.
\end{proof}
\begin{figure}
\begin{center}
\begin{tabular}{ccc}
 \scalebox{2.5}{\includegraphics{./tree-eps-converted-to.pdf}}
 &\qquad&
 \scalebox{2.5}{\includegraphics{./loop+tree-eps-converted-to.pdf}}
\end{tabular}
\end{center}
\caption{Diagrams without loop (left) and with a loop (right)}
\label{fig:loop+tree}
\end{figure}

The last column in Table \ref{tab:derivatives} shows the number of loops in the diagram.

\begin{table
\caption{Derivatives of correlation functions $\kappa^{(N)}$, $N=2,3,4$.}
\label{tab:derivatives}
\begin{center}
\begin{tabular}{ccccc}
\hline
Abbr. & \parbox[c][6ex][c]{12ex}{Representative \\ derivative} & Multiplicity & Diagram & \parbox[c][6ex][c]{5ex}{\# of \\ loops} \\ \hline
 $\rho_1$     & $\rho'(0)$ & 1 &
 \scalebox{1.5}{\includegraphics{./k1-eps-converted-to.pdf}} & 0 \\
 $\rho_2$     & $\rho''(0)$ & 1 &
 \scalebox{1.5}{\includegraphics{./k2-eps-converted-to.pdf}} & 1 \\
\hline
 $\kappa_0$     & $\kappa^{(3)}(0,0,0)$ & 1 & & 0 \\
 $\kappa_1$    & $\kappa^{(3)}_{(1,0,0)}(0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./k1-eps-converted-to.pdf}} & 0 \\
 $\kappa_2$    & $\kappa^{(3)}_{(2,0,0)}(0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./k2-eps-converted-to.pdf}} & 1 \\
 $\kappa_{11}$ & $\kappa^{(3)}_{(1,1,0)}(0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./k11-eps-converted-to.pdf}} & 0 \\
 $\kappa_{21}$ & $\kappa^{(3)}_{(2,1,0)}(0,0,0)$ & 6 &
 \scalebox{1.5}{\includegraphics{./k21-eps-converted-to.pdf}} & 1 \\
 $\kappa_{111}$ & $\kappa^{(3)}_{(1,1,1)}(0,0,0)$ & 1 &
 \scalebox{1.5}{\includegraphics{./k111-eps-converted-to.pdf}} & 1 \\
\hline
 $\kk_0$        & $\kappa^{(4)}_{(0,0,0,0,0,0)}(0,0,0,0,0,0)$ & 1 &
 & 0 \\
 $\kk_1$        & $\kappa^{(4)}_{(1,0,0,0,0,0)}(0,0,0,0,0,0)$ & 6 &
 \scalebox{1.5}{\includegraphics{./k1-eps-converted-to.pdf}} & 0 \\
 $\kk_2$        & $\kappa^{(4)}_{(2,0,0,0,0,0)}(0,0,0,0,0,0)$ & 6 &
 \scalebox{1.5}{\includegraphics{./k2-eps-converted-to.pdf}} & 1 \\
 $\kk_{11}^a$   & $\kappa^{(4)}_{(1,1,0,0,0,0)}(0,0,0,0,0,0)$ & 12 &
 \scalebox{1.5}{\includegraphics{./k11-eps-converted-to.pdf}} & 0 \\
 $\kk_{11}^{aa}$& $\kappa^{(4)}_{(1,0,0,0,0,1)}(0,0,0,0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./l11aa-eps-converted-to.pdf}} & 0 \\
 $\kk_{111}^c$  & $\kappa^{(4)}_{(1,1,0,1,0,0)}(0,0,0,0,0,0)$ & 4 &
 \scalebox{1.5}{\includegraphics{./k111-eps-converted-to.pdf}} & 1 \\
 $\kk_{111}^d$  & $\kappa^{(4)}_{(1,1,1,0,0,0)}(0,0,0,0,0,0)$ & 4 &
 \scalebox{1.5}{\includegraphics{./l111d-eps-converted-to.pdf}} & 0 \\
 $\kk_{111}^a$  & $\kappa^{(4)}_{(1,1,0,0,1,0)}(0,0,0,0,0,0)$ & 12 &
 \scalebox{1.5}{\includegraphics{./l111a-eps-converted-to.pdf}} & 0 \\
 $\kk_{21}^b$   & $\kappa^{(4)}_{(2,1,0,0,0,0)}(0,0,0,0,0,0)$ & 24 &
 \scalebox{1.5}{\includegraphics{./k21-eps-converted-to.pdf}} & 1 \\
 $\kk_{21}^{ab}$& $\kappa^{(4)}_{(2,0,0,0,0,1)}(0,0,0,0,0,0)$ & 6 &
 \scalebox{1.5}{\includegraphics{./l21ab-eps-converted-to.pdf}} & 1 \\
 $\kk_{22}^{bb}$& $\kappa^{(4)}_{(2,0,0,0,0,2)}(0,0,0,0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./l22bb-eps-converted-to.pdf}} & 2 \\
 $\kk_{211}^f$  & $\kappa^{(4)}_{(2,1,0,0,1,0)}(0,0,0,0,0,0)$ & 12 &
 \scalebox{1.5}{\includegraphics{./l211f-eps-converted-to.pdf}} & 1 \\
 $\kk_{211}^{b}$& $\kappa^{(4)}_{(2,1,0,0,0,1)}(0,0,0,0,0,0)$ & 24 &
 \scalebox{1.5}{\includegraphics{./l211b-eps-converted-to.pdf}} & 1 \\
 $\kk_{211}^d$  & $\kappa^{(4)}_{(2,1,1,0,0,0)}(0,0,0,0,0,0)$ & 12 &
 \scalebox{1.5}{\includegraphics{./l211d-eps-converted-to.pdf}} & 1 \\
 $\kk_{1111}^c$ & $\kappa^{(4)}_{(1,1,0,0,1,1)}(0,0,0,0,0,0)$ & 3 &
 \scalebox{1.5}{\includegraphics{./l1111c-eps-converted-to.pdf}} & 1 \\
 $\kk_{1111}^q$ & $\kappa^{(4)}_{(1,1,1,1,0,0)}(0,0,0,0,0,0)$ & 12 &
 \scalebox{1.5}{\includegraphics{./l1111q-eps-converted-to.pdf}} & 1 \\
\hline
\end{tabular}
\end{center}
\end{table}

Using the abbreviated notations listed in Table \ref{tab:derivatives},
the cumulants $K_{u,v,w}(s,T,\Theta)$ in (\ref{Kuvw})
 are evaluated.
As an example, we consider
\begin{align*}
\nu K_{0,2,1}(s,T,\Theta)
=& \mathrm{cum}(\langle T,\nabla X\rangle,\langle T,\nabla X\rangle,\mathrm{tr}(\Theta \nabla^2 X))\\
=& \sum_{i,j,k,l=1}^n \tau_i \tau_j \theta_{ij} \mathrm{cum}(X_i,X_j,X_{kl}),
\end{align*}
and
\begin{align*}
& \mathrm{cum}(X_i,X_j,X_{kl}) \\
&= \frac{\partial^4}{\partial t_1^i\partial t_2^j\partial t_3^k\partial t_3^l}
\nu \kappa^{(3)}\bigl(\tfrac12\Vert t_1-t_2\Vert^2,\tfrac12\Vert t_1-t_3\Vert^2,\tfrac12\Vert t_2-t_3\Vert^2\bigr)
 \big|_{t_1=t_2=t_3} \\
&= \nu \bigl(-2\kappa_{11} \delta_{ij}\delta_{kl} + \kappa_{11} (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\bigr),
\end{align*}
hence we obtain the result
$K_{0,2,1}(s,T,\Theta)
= -2\kappa_{11}\Vert T\Vert^2\mathrm{tr}(\Theta) + 2 \kappa_{11} T^\top\Theta T$.
This manipulation can be conducted by the symbolic computation.
We list the cumulants $K_{u,v,w}(s,T,\Theta)$ up to the 4-th order below.

\subsubsection*{The 2nd cumulants.}
\begin{align*}
& K_{2,0,0}
 = s^2, \qquad
K_{0,2,0}
 = \gamma \Vert T\Vert^2, \quad
K_{1,1,0} = 0, \\
& \mathrm{cum}(\mathrm{tr}(\Theta R),\mathrm{tr}(\Theta R))
= \alpha \mathrm{tr}(\Theta^2) + \beta \mathrm{tr}(\Theta)^2, \\
& \mathrm{cum}(s X,\mathrm{tr}(\Theta R)) =
 \mathrm{cum}(\langle T,\nabla X\rangle,\mathrm{tr}(\Theta R)) = 0,
\end{align*}
where
$\alpha = 2\rho_2$, and
$\beta = \rho_2-\rho_1^2$.

\subsubsection*{The 3rd cumulants:}
\begin{align*}
 K_{3,0,0}
 =& C_0 s^3, \quad
K_{1,2,0}
 = C_1 s \Vert T\Vert^2, \quad
K_{2,0,1}
= s^2 F_0 \mathrm{tr}(\Theta), \\
K_{1,0,2}
 =& F_1 s \mathrm{tr}(\Theta)^2 + F_2 s \mathrm{tr}(\Theta^2), \quad
K_{0,2,1}
 = J_1 \Vert T\Vert^2 \mathrm{tr}(\Theta) + J_2 T^\top\Theta T, \\
K_{1,1,1}
 =& K_1 \mathrm{tr}(\Theta)^3 + K_2 \mathrm{tr}(\Theta^2)\mathrm{tr}(\Theta) + K_3 \mathrm{tr}(\Theta^3),
\end{align*}
where
\begin{align*}
& C_0 = \kappa_0, \quad
  C_1 = -\kappa_1, \quad
  F_0 = 2\kappa_1, \quad
  F_1 = 3\kappa_{11}+\kappa_2, \quad
  F_2 = 2\kappa_2, \\
& J_1 = -2\kappa_{11}, \quad
  J_2 = 2\kappa_{11}, \quad
 K_1 = 2(\kappa_{111}+3\kappa_{21}), \quad
 K_2= 12\kappa_{21}, \quad
 K_3= -8\kappa_{111}.
\end{align*}

\subsubsection*{The 4th cumulants.}
\begin{align*}
K_{4,0,0}
 =& \CC_0 s^4, \quad
K_{2,2,0}
 = \CC_1 s^2 \Vert T\Vert^2, \quad
K_{0,4,0}
 = \CC_2 \Vert T\Vert^4, \quad
K_{3,0,1}
 = \FF_0 s^3 \mathrm{tr}(\Theta), \\
K_{2,0,2}
 =& \FF_1 s^2 \mathrm{tr}(\Theta)^2 + \FF_2 s^2 \mathrm{tr}(\Theta^2), \quad
K_{1,2,1}
 = \JJ_1 s \Vert T\Vert^2 \mathrm{tr}(\Theta) + \JJ_2 s T^\top\Theta T, \\
K_{1,0,3}
 =& \KK_1 s \mathrm{tr}(\Theta)^3 + \KK_2 s \mathrm{tr}(\Theta^2)\mathrm{tr}(\Theta) + \KK_3 s \mathrm{tr}(\Theta^3), \\
K_{0,2,2}
 =& \LL_1 \Vert T\Vert^2 \mathrm{tr}(\Theta)^2 + \LL_2 T^\top\Theta T \mathrm{tr}(\Theta) + \LL_3 \Vert T\Vert^2 \mathrm{tr}(\Theta^2) + \LL_4 T^\top\Theta^2 T, \\
K_{0,0,4}
 =& \MM_1 \mathrm{tr}(\Theta)^4 + \MM_2 \mathrm{tr}(\Theta^2)\mathrm{tr}(\Theta)^2 + \MM_3 \mathrm{tr}(\Theta^2)^2 + \MM_4 \mathrm{tr}(\Theta^3)\mathrm{tr}(\Theta) + \MM_5 \mathrm{tr}(\Theta^4),
\end{align*}
where
\begin{align*}
& \CC_0 = \kk_0, \quad
  \CC_1 = -\kk_1, \quad
  \CC_2 = 3 \kk_{11}^{aa}, \quad
  \FF_0 = 3\kk_1, \\
& \FF_1 = 6\kk_{11}^a+2\kk_{11}^{aa}+\kk_2, \quad
  \FF_2 = 2\kk_2, \quad
  \JJ_1 = -(2\kk_{11}^a +\kk_{11}^{aa}), \quad
  \JJ_2 = 2\kk_{11}^a, \\
& \KK_1 = 12\kk_{111}^a+2\kk_{111}^c+4\kk_{111}^d+3\kk_{21}^{ab}+6\kk_{21}^b, \quad
  \KK_2 = 6(\kk_{21}^{ab} + 2\kk_{21}^b), \\
& \KK_3 = -8\kk_{111}^c, \quad
  \LL_1 = -(6\kk_{111}^a +2\kk_{111}^d +\kk_{21}^{ab}), \quad
  \LL_2 = 4(2\kk_{111}^a +\kk_{111}^d), \\
& \LL_3 = -2\kk_{21}^{ab}, \quad
 \LL_4 = -8\kk_{111}^a, \\
& \MM_1 = 3(2\kk_{1111}^c +8\kk_{1111}^q +8\kk_{211}^b +4\kk_{211}^d +4\kk_{211}^f +\kk_{22}^{bb}) \\
& \MM_2 = 12(4\kk_{211}^b+2\kk_{211}^d+2\kk_{211}^f+\kk_{22}^{bb}), \quad
  \MM_3 = 12\kk_{22}^{bb}, \quad
  \MM_4 = -96\kk_{1111}^q, \\
& \MM_5 = 48\kk_{1111}^c.
\end{align*}

\subsection{Perturbation of the Euler characteristic density}
\label{sec:EC_density}

This section attempts to obtain the perturbation formula for the Euler characteristic density in (\ref{Xin}), which is expressed by $(X,\nabla X,R)$.
In the previous section, we obtained the Taylor series of the cumulant generating function $\log\psi(s,T,\Theta)$ in (\ref{mgf_X''}) up to $O(\nu^2)$.
The moment generating function of $(X,\nabla X,R)$ was $\psi(s+\gamma\mathrm{tr}(\Theta),T,\Theta)$.
We will conduct the formal calculation, however, its validity condition is provided in Section \ref{sec:regularity}.

The cumulant generating function is
\begin{align*}
\log\psi(\s,T,\Theta)
=& \log\psi^0_{X}(s) + \log\psi^0_{V}(T) +\log\psi^0_{R}(\Theta) \\
 & + \nu Q_{3,n}(\s,T,\Theta) + \nu^2 Q_{4,n}(\s,T,\Theta) + o(\nu^2),
\end{align*}
where $\s=\s(s,\Theta)=s+\gamma\mathrm{tr}(\Theta)$,
\begin{equation*}
 \psi^0_{X}(s) = e^{\frac{1}{2}s^2}, \qquad
 \psi^0_{V}(T) = e^{\frac{1}{2}\gamma \Vert T\Vert^2}, \qquad
 \psi^0_{R}(\Theta) = e^{\frac{1}{2}\alpha\mathrm{tr}(\Theta^2) +\frac{1}{2}\beta\mathrm{tr}(\Theta)^2},
\end{equation*}
and
\begin{align*}
 Q_{k,n}(\s,T,\Theta) = \sum_{u+v+w=k}\frac{1}{u!\,v!\,w!} K_{u,v,w}(\s,T,\Theta),\ \ k=3,4.
\end{align*}
Therefore,
\begin{align}
\label{mgf}
 \psi(\s,T,\Theta)
= \psi_X^0(s)\psi_V^0(T)\psi_R^0(\Theta)
 \Bigl(1 +\nu Q_{3,n}(\s,T,\Theta) +\nu^2 Q_{4,n}(\s,T,\Theta) \\ +\frac{1}{2}\nu^2 Q_{3,n}(\s,T,\Theta)^2 + o(\nu^2) \Bigr). \nonumber
\end{align}
In the following calculation, we will ignore the term $o(\nu^2)$.
This formal manipulation is validated in Section \ref{sec:regularity}.

Define
\begin{align}
\label{psi_V0}
 \psi_{V=0}(\s,\Theta)
=& \int_{\mathbb{R}} \int_{\mathrm{Sym}(n)} e^{s x+\mathrm{tr}(\Theta R)} p(x,0,R) \dd R \dd x \\
=& \frac{1}{(2\pi)^n}\int_{\mathbb{R}^n} \psi(\s,i T,\Theta) \dd T. \nonumber
\end{align}
Recall that the general term of the expansion in (\ref{mgf}) is (\ref{K}).
The integration with respect to $\dd T$ is conducted as an expectation with respect to the Gaussian random variable $T\sim N_n(0,\gamma^{-1})$,
and multiplying it by the normalizing factor
\begin{equation}
\label{normatizing-T}
 \frac{1}{(2\pi)^n}\int_{\mathbb{R}^n}e^{-\frac{1}{2}\gamma\Vert T\Vert^2} \dd T = (2\pi\gamma)^{-n/2}.
\end{equation}
Let
\[
 (n)_m = n(n-1)\cdots(n-m+1) = \frac{\Gamma(n+1)}{\Gamma(n-m+1)}
\]
be the falling factorial.
Since $\Vert T\Vert$ is independent of
$\prod_{j\ge 1} (T^\top\Theta^j T/\Vert T\Vert^2)^{v_j}$,
the expectation of the factor of (\ref{K}) including $T$ becomes
\begin{align}
\label{ET}
 \mathbb{E}^{T}\biggl[\prod_{j\ge 0}(T^\top \Theta^j T)^{v_j}\biggr]
=& \mathbb{E}^{T}\bigl[\Vert T\Vert^{2 v}\bigr] \mathbb{E}^{T}\biggl[\frac{\prod_{j\ge 1}(T^\top \Theta^j T)^{v_j}}{\Vert T\Vert^{2(v-v_0)}}\biggr] \\
=& \mathbb{E}^{T}\bigl[\Vert T\Vert^{2 v}\bigr] \frac{\mathbb{E}^{T}\bigl[\prod_{j\ge 1}(T^\top \Theta^j T)^{v_j}\bigr]}{\mathbb{E}^{T}\bigl[\Vert T\Vert^{2(v-v_0)}\bigr]} \nonumber \\
=& (-2/\gamma)^{v_0} (-(n/2+v-v_0))_{v_0} \times \zeta_{v_1,v_2,\ldots}(\Theta), \nonumber
\end{align}
where $v=\sum_{j\ge 0} v_j$, and
\begin{equation*}
 \zeta_{v_1,v_2,\ldots}(\Theta)
 = \mathbb{E}\biggl[\prod_{j\ge 1}(\xi^\top \Theta^j \xi)^{v_j}\biggr], \quad \xi\sim N_n(0,I_n)
\end{equation*}
is a polynomial in $\mathrm{tr}(\Theta^j)$, $j\ge 1$.
For example, $\zeta_{1}(\Theta)=\mathrm{tr}(\Theta)$, $\zeta_{2}(\Theta)=\mathrm{tr}(\Theta^2)$, $\zeta_{1,1}(\Theta)=\mathrm{tr}(\Theta)^2+2\mathrm{tr}(\Theta^2)$.
Note that $\zeta_{v_1,v_2,\ldots}(\Theta)$ does not include the dimension $n$ explicitly.

For the terms in $Q_3$, $Q_4$ and $\frac{1}{2}Q_3^2$ in (\ref{mgf}) including $T$,
we can rewrite the terms according to the rules below:
\begin{align}
\begin{aligned}
 \Vert T\Vert^{2} &\mapsto \gamma^{-1} n, \qquad
 \Vert T\Vert^{4} \mapsto \gamma^{-2} n(n+2), \\
 T^\top\Theta^k T &\mapsto \gamma^{-1} \mathrm{tr}(\Theta^k), \quad k=1,2, \\
 \Vert T\Vert^2 \cdot T^\top\Theta T &\mapsto \gamma^{-2}(n+2)\mathrm{tr}(\Theta), \\
 (T^\top\Theta T)^2 &\mapsto \gamma^{-2}[\mathrm{tr}(\Theta)^2+2\mathrm{tr}(\Theta^2)],
\end{aligned}
\label{rule1}
\end{align}
and multiply the normalizing factor (\ref{normatizing-T}).
Then, we have
\begin{align}
\label{mgf_V}
 \psi_{V=0}(\s,\Theta)
 = \psi^0_{X}(s) \psi^0_{R}(\Theta) (2\pi\gamma)^{-n/2}
  \Bigl(1 + \nu\widetilde Q_{3,n}(\s,\Theta) + \nu^2\widetilde Q_{4,n}(\s,\Theta) \\
 + \frac{1}{2}\nu^2\widetilde Q_{3,n}^{(2)}(\s,\Theta)
 + o(\nu^2) \Bigr), \nonumber
\end{align}
where
$\widetilde Q_{3,n}(\s,\Theta)$, $\widetilde Q_{4,n}(\s,\Theta)$, and $\widetilde Q_{3,n}^{(2)}(\s,\Theta)$
are
$Q_{3,n}(\s,T,\Theta)$, $Q_{4,n}(\s,T,\Theta)$, $Q_{3,n}(\s,T,\Theta)^2$ with terms including $T$ replaced according to the rules (\ref{rule1}).

Recall that the integral we are in the process of obtaining is $\Xi_n(v)$ in (\ref{Xin}).
From the definition of $\psi_{V=0}(i \s,\Theta)$ in (\ref{psi_V0}),
we see that (\ref{Xin}) is evaluated as
\begin{equation}
\label{target}
 \Xi_n(x) = \int_{x}^\infty \biggl[
 \det(-D_\Theta + \gamma x I_n) \Big|_{\Theta=0}
 \frac{1}{2\pi}\int_{\mathbb{R}} e^{-isx} \psi_{V=0}(i \s,\Theta) \dd s \biggr] \dd x,
\end{equation}
where $D_\Theta$ is an $n\times n$ symmetric matrix differential operator defined by
\begin{equation}
\label{DTheta}
(D_\Theta)_{ij} = \frac{1+\delta_{ij}}{2}\frac{\partial}{\partial (\Theta)_{ij}} = \frac{\partial}{\partial \tau_{ij}} \quad (i\le j).
\end{equation}

The integral in (\ref{target}) with respect to $\dd s$ is easily evaluated by
\begin{equation}
 \frac{1}{2\pi}\int_{\mathbb{R}} e^{-isx} \psi^0_{X}(s) (is)^k \dd s = H_k(x) \phi(x),
\label{s2x}
\end{equation}
where
\[
 \phi(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}
\]
is the probability density of the standard Gaussian distribution $N(0,1)$, and
\[
 H_k(x) = \phi(x)^{-1}(-\dd/\dd x)^k\phi(x)
\]
is the Hermite polynomial of degree $k$.
For the derivatives with respect to $\Theta$, we use the lemma below.
The main part of the proof is given in Section \ref{sec:Hermite}.

\begin{proposition}
\label{prop:DThetadet}
\begin{align}
\label{theta2x}
 \det(-D_\Theta + \gamma x I_n) \bigl(\psi_R^0(\Theta)\mathrm{tr}(\Theta^{c_1})\cdots\mathrm{tr}(\Theta^{c_k})\bigr)\Big|_{\Theta=0} \\
= (-1)^m \gamma^{n-m} (-1/2)^{m-k} (n)_m H_{n-m}(x), \nonumber
\end{align}
where $m=\sum_{i=1}^k c_i$.
\end{proposition}

\begin{proof}
Recall that $\gamma=\sqrt{\alpha/2-\beta}$.
(\ref{theta2x2}) in Lemma \ref{lem:DThetadet2} with $\beta:=\beta/\gamma^2$, $\Theta:=\gamma\Theta$, and $x:=-x$ yields (\ref{theta2x}).
\end{proof}

Now, we are ready to complete the whole calculations.
\begin{itemize}
\item[Step 1.]
Expand the inside of the parenthesis in (\ref{mgf_V}).
The resulting mgf will be $\psi_X^0(s)\psi_R^0(\Theta)$ times polynomial in $s$ and $\mathrm{tr}(\Theta^k)$, $k\le 4$.

\item[Step 2.]
Applying (\ref{s2x}) and (\ref{theta2x}) yields the result of Step 1,
which is $\phi(x)$ multiplied by a polynomial in $x^k$, and $H_{k}(x)$.

\item[Step 3.]
By using the three-term relation
\begin{equation*}
 x H_{k+1}(x) = H_{k+1}(x) + k H_{k-1}(x),
\end{equation*}
we arrange the result of Step 2 to be $\phi(x)$ multiplied by a linear combination of $H_{k}(x)$.

\item[Step 4.]
Taking the integration $\int_{x}^\infty \dd x$ by using
\[
 \int_{x}^\infty H_k (x) \phi(x) \dd x = H_{k-1}(x) \phi(x),
\]
we complete the integral (\ref{target}).
Let $H_{-1}(x) = \phi(x)^{-1}\int_x^\infty \phi(x) \dd x$.
\end{itemize}
These processes are simple rewritings of terms that can be conducted by symbolic calculation.
The results are summarized in the theorem below.

\begin{theorem}
\label{thm:perturbation}
Under Assumptions \ref{as:B-R} and \ref{as:A-T}, it holds that
\begin{align}
\label{Xi_perturbation}
 \gamma^{-n/2}\Xi_n(x)
=& (2\pi)^{-n/2} \phi(x) \times\Bigl( H_{n-1}(x) + \nu \Delta_{1,n}(x) + \nu^2 \Delta_{2,n}(x) \Bigr) + o(\nu^2),
\end{align}
where
\begin{align}
\label{D1n}
\Delta_{1,n}(x) =
&  \tfrac{1}{2} \gamma^{-2} \kappa_{11} n (n-1) H_{n-2}(x)
 - \tfrac{1}{2} \gamma^{-1} \kappa_1 n H_{n}(x)
 + \tfrac{1}{6} \kappa_0 H_{n+2}(x),
\end{align}
\begin{align}
\label{D2n}
\Delta_{2,n}(x) =
& \Bigl( -\tfrac{1}{6}\gamma^{-3} (3\kk_{111}^a + \kk_{111}^d)
 +\tfrac{1}{8}\gamma^{-4} \kappa_{11}^2 (n-7) \Bigr) n(n-1)(n-2) H_{n-3}(x) \nonumber \\
& + \Bigl( \tfrac{1}{8}\gamma^{-2} \bigl(\kk_{11}^{aa} (n-2) +4\kk_{11}^a (n-1)\bigr) \nonumber \\
&\quad -\tfrac{1}{4}\gamma^{-3} \kappa_{1}\kappa_{11} (n-1)(n-4) \Bigr) n H_{n-1}(x) \nonumber \\
& + \Bigl( - \tfrac{1}{4}\gamma^{-1} \kk_{1}
 +\tfrac{1}{24}\gamma^{-2} \bigl(3\kappa_1^2 (n-2) + 2\kappa_0\kappa_{11} (n-1)\bigr) \Bigr) n H_{n+1}(x) \nonumber \\
& + \Bigl( \tfrac{1}{24} \kk_0 -\tfrac{1}{12}\gamma^{-1} \kappa_{0}\kappa_{1} n \Bigr) H_{n+3}(x)
 +\tfrac{1}{72} \kappa_{0}^2 H_{n+5}(x). \nonumber
\end{align}
\end{theorem}
\begin{proof}
For the regularity conditions, see Section \ref{sec:regularity}.
\end{proof}

This theorem shows the kind of non-Gaussianity captured by the Euler characteristic/Minkowski functional approach.
Formula (\ref{Xi_perturbation}) includes
\begin{quote}
$\rho_0=1$, $\rho_1=-\gamma$, $\kappa_0$, $\kappa_1$, $\kappa_{11}$,
$\kk_0$, $\kk_1$, $\kk_{11}^a$, $\kk_{11}^{aa}$, $\kk_{111}^a$, $\kk_{111}^d$,
\end{quote}
and excludes
\begin{quote}
$\rho_2$, $\kappa_2$, $\kappa_{21}$, $\kappa_{111}$,
$\kk_2$, $\kk_{111}^c$, $\kk_{21}^b$, $\kk_{21}^{ab}$, $\kk_{22}^{bb}$, $\kk_{211}^f$, $\kk_{211}^{b}$, $\kk_{211}^d$, $\kk_{1111}^c$, $\kk_{1111}^q$.
\end{quote}
According to Table \ref{tab:derivatives},
the diagrams of the excluded parameters contain loops, while those of the included parameters are loop-free.
Theorem \ref{thm:loop} claims that this observation holds for any cumulant $\kappa^{(N)}$ in the perturbation of an arbitrary order.
Let
\[
 \Pi_K(\Theta) = \mathrm{tr}(\Theta)^K - (-2)^{K-1}\mathrm{tr}(\Theta^K) \quad (K\ge 2).
\]

\begin{theorem}
\label{thm:loop}
Suppose that the diagram of the derivative of the $N$-point function $\kappa^{(N)}$,
\begin{equation}
\label{kappaN1}
 \biggl(\prod_{(a,b)}\frac{\partial}{\partial x_{ab}}\biggr)\kappa^{(N)}(0,\ldots,0),
\end{equation}
has a loop of length $K\ge 2$,
where $\prod_{(a,b)}$ runs over all edges of the diagram.
Then, (\ref{kappaN1}) does not appear in the asymptotic expansion of arbitrary order.
\end{theorem}

\begin{proof}
We first assume that the diagram is connected.
Let $\widetilde K \,(\le N)$ be the number of vertices of the diagram.
When there exists a loop in the diagram, as seen in Lemma \ref{lem:diagram},
the edge set of the diagram is disjointly divided as $L\sqcup T$, where
\[
 L = \{(1,2),(2,3),\ldots,(K-1,K),(K,1) \}
\]
is a loop of length $K$, and $T$ is a union of the tree edges.
For example, in Figure \ref{fig:loop+tree} (right), $K=3$, $\widetilde K=6$,
$L=\{(1,2),(2,3),(1,3)\}$, and $T=\{(1,4),(4,5),(4,6)\}$.
When a loop is included, as in the proof of Lemma \ref{lem:diagram},
the derivative (\ref{kappaN1}) appears only from the operation
\begin{equation}
\label{taylor0}
 \prod_{a=1}^{\widetilde K} \mathrm{tr}(\Theta\nabla_{t_a}^2)
 \kappa^{(N)}\bigl((\tfrac{1}{2}\Vert t_a-t_b\Vert^2)_{a<b}\bigr)\Big|_{t_1=\cdots=t_N}.
\end{equation}
By picking up the non-vanishing term in the Taylor expansion of $\kappa^{(N)}$,
(\ref{taylor0}) is proportional to
\begin{equation}
\label{taylor1}
 \prod_{a=1}^{\widetilde K} \mathrm{tr}(\Theta\nabla_{t_a}^2)
 \prod_{(a,b)\in L\sqcup T}\bigl(\tfrac{1}{2}\Vert t_a-t_b\Vert^2\bigr)\Big|_{t_1=\cdots=t_N}\times \biggl(\prod_{(a,b)\in L\sqcup T} \frac{\partial}{\partial x_{ab}}\biggr)\kappa^{(N)}(0,\ldots,0).
\end{equation}
We shall prove later that (\ref{taylor1}) is equal to
\begin{equation}
\label{coefficient}
 2 \mathrm{tr}(\Theta)^{\widetilde K-K}\Pi_{K}(\Theta)
 \biggl(\prod_{(a,b)\in L\sqcup T} \frac{\partial}{\partial x_{ab}}\biggr)\kappa^{(N)}(0,\ldots,0).
\end{equation}

When the diagram has two or more connecting components,
the coefficient in (\ref{coefficient}) is multiplied by the coefficient for the other components.
The multiplying factor is a polynomial in $T^\top\Theta^k T$ ($k\ge 0$) and $\mathrm{tr}(\Theta^k)$ ($k\ge 1$).
After taking the expectation with respect to $T\sim N_n(0,\gamma^{-1}I_n)$, they will be a polynomial in $\mathrm{tr}(\Theta^k)$ ($k\ge 1$).
Therefore, if the diagram has a loop of length $K$, after taking the expectation with respect to $T$, the coefficient of (\ref{kappaN1}) will be a polynomial in $\mathrm{tr}(\Theta^k)$ ($k\ge 1$) with a factor $\Pi_K(\Theta)$.
It will vanish at the final stage because of Proposition \ref{prop:DThetadet},
\begin{align*}
\det & (-D_\Theta + \gamma x I_n) \bigl(\psi_R^0(\Theta)\mathrm{tr}(\Theta^{c_1})\cdots\mathrm{tr}(\Theta^{c_k}) \Pi_K(\Theta) \bigr)\Big|_{\Theta=0} \\
=& (-1)^{m+K} \gamma^{n-(m+K)} (-1/2)^{(m+K)-(k+K)} (n)_{m+K} H_{n-(m+K)}(x) \\
& -(-2)^{K-1} (-1)^{m+K} \gamma^{n-(m+K)} (-1/2)^{(m+K)-(k+1)} (n)_{m+K} H_{n-(m+K)}(x) \\
=& 0,
\end{align*}
where $m=\sum_{i=1}^k c_i$.

Finally we prove (\ref{coefficient}).
First, since $T$ is a union of trees,
\begin{equation}
\label{LT}
 \biggl(\prod_{a=1}^{\widetilde K} \frac{\partial^2}{\partial t_{a}^{i_a}\partial t_{a}^{j_a}}\biggr)
 \Delta_L \Delta_T
 = \biggl(\prod_{a=1}^{K} \frac{\partial^2}{\partial t_{a}^{i_a}\partial t_{a}^{j_a}}\biggr) \Delta_L \prod_{a=K+1}^{\widetilde K} \delta_{i_a j_a},
\end{equation}
where
\[
 \Delta_L = \prod_{(a,b)\in L} \bigl(\tfrac{1}{2}\Vert t_a-t_b\Vert^2 \bigr), \qquad
 \Delta_T = \prod_{(a,b)\in T} \bigl(\tfrac{1}{2}\Vert t_a-t_b\Vert^2 \bigr).
\]
Let $u_{a}^{i} = t_{a}^{i}-t_{a-1}^{i}$ ($t_{0}^{i}=t_{K}^{i}$).
Since $\partial/\partial t_{a}^{i}=\partial/\partial u_{a}^{i} - \partial/\partial u_{a-1}^{i}$,
\begin{align}
\label{L}
\prod_{a=1}^K \frac{\partial^2}{\partial t_{a}^{i_a}\partial t_{a}^{j_a}}\Delta_L
&= \prod_{a=1}^K \biggl(\frac{\partial^2}{\partial u_{a}^{i}\partial u_{a}^{j}} + \frac{\partial^2}{\partial u_{a-1}^{i}\partial u_{a-1}^{j}}
   - \frac{\partial^2}{\partial u_{a}^{i}\partial u_{a-1}^{j}} - \frac{\partial^2}{\partial u_{a-1}^{i}\partial u_{a}^{j}}\biggr)\Delta_L \\
&= 2\biggl(\prod_{a=1}^K \frac{\partial^2}{\partial u_{a}^{i_a}\partial u_{a}^{j_a}} \biggr)\Delta_L
- 2\biggl(\sum_{\rm chain} \prod_{a=1}^K(-1)\frac{\partial^2}{\partial u_{a}^{i_a}\partial u_{a}^{j_a}}\biggr)\Delta_L, \nonumber
\end{align}
where $\sum_{\rm chain}$ runs over the set such that $\{(i_1,j_1),\ldots,(i_K,j_K)\}$ forms an undirected chain.
Hence, we have
\begin{align*}
\biggl(\prod_{a=1}^{\widetilde K} \mathrm{tr}(\Theta\nabla_{t_a}^2)\biggr) \Delta_L \Delta_T
=& \sum_{{i_1,\ldots,i_{\widetilde K},j_1,\ldots,j_{\widetilde K}}=1}^n
 \theta_{i_1 j_1}\cdots\theta_{i_{K'} j_{K'}}
 \biggl(\prod_{a=1}^{\widetilde K} \frac{\partial^2}{\partial t_{a}^{i_a}\partial t_{a}^{j_a}}\biggr)\Delta_L \Delta_T \\
=& 2\mathrm{tr}(\Theta)^{\widetilde K-K}\bigl(\mathrm{tr}(\Theta)^K - (-2)^K \mathrm{tr}(\Theta^K)\bigr)
\end{align*}
because of (\ref{LT}) and (\ref{L}).
\end{proof}

\subsection{Expected Minkowski functional on a bounded domain}
\label{sec:boundary}

In the previous subsection, we have obtained the perturbation formula for the Euler characteristic density $\Xi_n(v)$ in (\ref{Xin}).
Here, we evaluate $\Xi_{d,e}(v)$ in (\ref{Xide}), and $\mathbb{E}[\chi(E_v)]$ in (\ref{LXi}), in turn.
They are shown to be expressed by the function $\Xi_n(v)$.
The argument in this subsection is valid for an arbitrary order in which the asymptotic expansion is validated.

The joint mgf of $(X,V_1,V_2,R_{11})$ expressing $\Xi_{d,e}(v)$ is
\begin{equation*}
 \mathbb{E}\left[\exp\Bigl( s X + (T_1^\top,T_2^\top) \begin{pmatrix}V_1 \\ V_2 \end{pmatrix} + \mathrm{tr}(\Theta_{11} R_{11}) \Bigr) \right]
 = \psi\biggl(\s,\begin{pmatrix}T_1 \\ T_2\end{pmatrix},\begin{pmatrix}\Theta_{11} & 0 \\ 0 & 0 \end{pmatrix}\biggr),
\end{equation*}
where $\psi$ is defined by (\ref{mgf_X''}) and $\s=s+\gamma\mathrm{tr}(\Theta_{11})$.

Define
\begin{align*}
 \psi_{V_1=0}(s,T_2,\Theta_{11})
 = \int_{\mathbb{R}}\int_{\mathrm{Sym}(d)}\int_{\mathbb{R}^{n-d}} & e^{s x+\langle T_2,V_2\rangle+\mathrm{tr}(\Theta_{11}R_{11})} \\
  & \times p(x,(0,V_2),R_{11}) \dd V_2 \dd R_{11} \dd x.
\end{align*}
Then,
\begin{equation}
\label{mgf_V1}
 \psi_{V_1=0}(\s,T_2,\Theta_{11})
 = \frac{1}{(2\pi)^d}\int_{\mathbb{R}^{d}}\psi\biggl(\s,\begin{pmatrix}i T_1 \\ T_2\end{pmatrix},\begin{pmatrix}\Theta_{11} & 0 \\ 0 & 0 \end{pmatrix}\biggr) \dd T_1.
\end{equation}
Using this, $\Xi_{d,e}(v)$ in (\ref{Xide}) is expressed as
\begin{equation}
\label{target1}
\Xi_{d,e}(v)
= \int_{x=v}^\infty \biggl[
\det\bigl(-D_{\Theta_{11}} + \gamma x I_{d-e}\bigr)\Big|_{\Theta_{11}=0} \frac{1}{2\pi}\int_{\mathbb{R}} e^{-i s x} \widetilde \psi_e(i \s,\Theta_{11}) \dd s \biggr] \dd x,
\end{equation}
where
\begin{align}
\widetilde \psi_e(\s,\Theta_{11})
= \int_{\mathbb{R}^{n-d}} \Vert V_2\Vert^e \frac{1}{(2\pi)^{n-d}}\int_{\mathbb{R}^{n-d}} e^{-i\langle T_2,V_2\rangle}\psi_{V_1=0}(\s,i T_{2},\Theta_{11}) \dd T_2 \dd V_{2}.
\label{psi_e}
\end{align}

In the following, we will establish the statement below.
\begin{lemma}
\label{lem:claim}
Under Assumption \ref{as:B-R'},
\begin{equation}
\label{claim}
\widetilde\psi_e(\s,\Theta_{11}) = \frac{\Omega_{n-d}}{\Omega_{n-d+e}} \times (\psi_{V=0}(\s,\Theta) \mbox{ with $n:=d-e$ and $\Theta:=\Theta_{11}$}) +o(\nu^{s_0-2}),
\end{equation}
where $\psi_{V=0}(\s,\Theta)$ is defined in (\ref{psi_V0}).
\end{lemma}

\begin{proof}
Under this assumption, the formal asymptotic expansion is validated up to $o(\nu^{s_0-2})$.
We evaluate (\ref{mgf_V1}) and (\ref{psi_e}).
There are three integrals with respect to $\dd T_1$, $\dd T_2$, and $\dd V_2$.
After expanding the $o(1)$ terms in the exponential function,
the integral (\ref{mgf_V1}) with respect to $\dd T_1$ becomes the expectation with respect to $T_1\sim N_{d}(0,\gamma^{-1}I_d)$,
and multiplying it by the factor
\[
 \frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^d} e^{-\frac{1}{2}\gamma\Vert T_1\Vert^2} \dd T_1 = (2\pi\gamma)^{-d/2}.
\]

For the integrals (\ref{psi_e}) with respect to $\dd T_2$ and $\dd V_2$, we need to evaluate
\[
 J_{e,k}(\gamma) = \int_{\mathbb{R}^{n-d}} \Vert V_2\Vert^e \frac{1}{(2\pi)^{n-d}}\int_{\mathbb{R}^{n-d}} e^{-i\langle T_2,V_2\rangle} e^{-\frac{1}{2}\gamma\Vert T_2\Vert^2} \Vert T_2\Vert^{2k} \dd T_2 \dd V_{2}.
\]
Noting that
\[
 \frac{1}{(2\pi)^{n-d}}\int_{\mathbb{R}^{n-d}} e^{-i\langle T_2,V_2\rangle} e^{-\frac{1}{2}\gamma\Vert T_2\Vert^2} \dd T_2 = \frac{1}{(2\pi\gamma)^{(n-d)/2}} e^{-\frac{1}{2\gamma}\Vert V_2\Vert^2},
\]
and
\begin{align*}
 J_{e,0}(\gamma)
 =& \int_{\mathbb{R}^{n-d}} \Vert V_2\Vert^e \frac{1}{(2\pi\gamma)^{(n-d)/2}} e^{-\frac{1}{2\gamma}\Vert V_2\Vert^2} \dd V_2 \\
=& \gamma^{e/2} \frac{2^{e/2}\Gamma(\frac{n-d+e}{2})}{\Gamma(\frac{n-d}{2})}
= (2\pi\gamma)^{e/2} \frac{\Omega_{n-d}}{\Omega_{n-d+e}},
\end{align*}
we have
\begin{align*}
 J_{e,k}(\gamma)/J_{e,0}(\gamma)
=& \Bigl(-2\frac{\partial}{\partial\gamma}\Bigr)^k J_{e,0}(\gamma)/J_{e,0}(\gamma) \\
=& \Bigl(-2\frac{\partial}{\partial\gamma}\Bigr)^k\gamma^{e/2}/\gamma^{e/2}
= (-2/\gamma)^k (e/2)_k.
\end{align*}
Therefore, rewriting the terms
\begin{equation}
 \Vert T_2\Vert^{2k} \mapsto (-2/\gamma)^k (e/2)_k, \quad
\label{rule3}
\end{equation}
and multiplying them by the factor $J_{e,0}(\gamma)$ yields the integral $J_{e,k}(\gamma)$.

When we consider the asymptotic expansion of an arbitrary order,
the general term has of the form of (\ref{K}), where $T=\begin{pmatrix} T_1 \\ T_2 \end{pmatrix}$, and $\Theta=\begin{pmatrix} \Theta_{11} & 0 \\ 0 & 0 \end{pmatrix}$.
The factor including $T_1$ and $T_2$ is
\[
 \prod_{j\ge 0} (T^\top\Theta^{j} T)^{v_j}
  = (\Vert T_1\Vert^2+\Vert T_2\Vert^2)^{v_0} \prod_{j\ge 1} (T_1^\top\Theta_{11}^j T_1)^{v_j}.
\]
We take the expectation under $T_1\sim N_d(0,\gamma^{-1}I_d)$.
As in the derivation of (\ref{ET}), the expectation of the right-hand side is
\begin{align*}
& \sum_{i=0}^{v_0} \binom{v_0}{i} \mathbb{E}^{T_1}\biggl[ \Vert T_1\Vert^{2(v_0-i)}
 \prod_{j\ge 1}(T_1^\top\Theta_{11}^{j}T_1)^{v_j} \biggr] \Vert T_2\Vert^{2 i} \\
&= \sum_{i=0}^{v_0} \binom{v_0}{i} (-2/\gamma)^{v_0-i}(-(d/2+v-v_0))_{v_0-i}\times \zeta_{v_1,v_2,\ldots}(\Theta_{11})\times \Vert T_2\Vert^{2 i},
\end{align*}
where $v=\sum_{j\ge 0}v_j$.
Applying the term-rewriting rule (\ref{rule3}) for $\Vert T_2\Vert^{2i}$, this becomes
\begin{align}
\label{result_de}
\sum_{i=0}^{v_0} \binom{v_0}{i} (-2/\gamma)^{v_0-i}\times (-(d/2+v-v_0))_{v_0-i} (-2/\gamma)^i (e/2)_i \times \zeta_{v_1,v_2,\ldots}(\Theta_{11}) \\
= (-2/\gamma)^{v_0} (-((d-e)/2+v-v_0))_{v_0}\times \zeta_{v_1,v_2,\ldots}(\Theta_{11}). \nonumber
\end{align}
Here we used the binomial formula for the falling factorial
\[
 (x+y)_{m} = \sum_{i=0}^m \binom{m}{i} (x)_{m-i} (y)_i.
\]
(\ref{result_de}) is equal to (\ref{ET}) with $n$ and $\Theta$ replaced by $d-e$ and $\Theta_{11}$, respectively.

Therefore, we have
\begin{align*}
& \widetilde\psi_e(\s,\Theta_{11})
= \psi^0_{X}(s) \psi^0_{R}(\Theta_{11}) (2\pi\gamma)^{-d/2} J_{e,0}(\gamma) \\
&\qquad\times \Bigl( 1 + \nu\widetilde Q_{3,d-e}(\s,\Theta_{11}) + \nu^2\widetilde Q_{4,d-e}(\s,\Theta_{11}) + \frac{1}{2}\nu^2\widetilde Q_{3,d-e}^{(2)}(\s,\Theta_{11}) + \cdots \Bigr),
\end{align*}
where $\widetilde Q_{3}$, $\widetilde Q_{4}$, and $\widetilde Q_{3}^{(2)}$ are defined in (\ref{mgf_V}).
Except for the factor $(2\pi\gamma)^{-e/2} J_{e,0}(\gamma) \allowbreak =\Omega_{n-d}/\Omega_{n-d+e}$,
$\widetilde\psi_e(\s,\Theta_{11})$ is exactly the same as the $\psi_{V=0}(\s,\Theta)$ in (\ref{mgf_V}) with $n:=d-e$ and $\Theta:=\Theta_{11}$.
\end{proof}

Therefore, from (\ref{claim}) combined with (\ref{target}) and (\ref{target1}),
we have
\[
 \Xi_{d,e}(v) = \frac{\Omega_{n-d}}{\Omega_{n-d+e}} \Xi_{d-e}(v) + o(\nu^{s_0-2}).
\]
Substituting this into (\ref{LXi}),
\begin{align*}
 \mathbb{E}[\chi(E_v)]
=& \sum_{d=0}^n \sum_{e=0}^d L_{n-d,e} \frac{\Omega_{n-d}}{\Omega_{n-d+e}} \Xi_{d-e}(v) + o(\nu^{s_0-2}) \\
=& \sum_{d=0}^n \Biggl(\sum_{e=d}^{n} L_{n-e,e-d} \frac{\Omega_{n-e}}{\Omega_{n-d}} \Biggr) \Xi_{d}(v) + o(\nu^{s_0-2}) \\
=& \sum_{d=0}^n \mathcal{L}_d \Xi_{d}(v) + o(\nu^{s_0-2}).
\end{align*}

\begin{theorem}
\label{thm:boundary}
Assume Assumptions \ref{as:B-R'} and \ref{as:A-T}.
Then,
\begin{equation}
\label{boundary}
 \mathbb{E}[\chi(E_v)] = \sum_{d=0}^n \mathcal{L}_d(E) \Xi_{d}(v) + o(\nu^{s_0-2}).
\end{equation}
$\mathcal{L}_d$ is the Lipschitz--Killing curvature of $E$ defined in (\ref{steiner}).
\end{theorem}

Note that if $p(x,V,R)$ has the moments of arbitrary order, then the remainder term of (\ref{boundary}) is $o(\nu^s)$, where $s$ is an arbitrary integer.

The Minkowski functional for the excursion set $\mathcal{M}_k(E_v)$ can be evaluated via the Crofton formula (Proposition \ref{prop:crofton}).
Let $L\in A(n,n-k)$, and let $X|_L$ be the restriction of $X$ on $L$, i.e., a random field on $E\cap L$.
Since $X$ is isotropic, when $L$ is given, we have
\[
  \mathbb{E}[\chi((E\cap L)_v)]
 = \sum_{d=0}^{n-k} \mathcal{L}_d(E\cap L)\,\Xi_{d}(v) + o(\nu^{s_0-2})\mathds{1}(E\cap L\ne\emptyset).
\]
We perform an integration with respect to $\int_{A(n,n-k)}d\mu_{n,n-k}(L)$.
Noting that $(E\cap L)_v=E_v\cap L$, we have
\[
  c_{n,0,k} \mathbb{E}[\mathcal{L}_{k}(E_v)]
 = \sum_{d=0}^{n-k} c_{n,d,k} \mathcal{L}_{k+d}(E)\,\Xi_{d}(v) + o(\nu^{s_0-2}) \omega_k.
\]
Let
\[
  \left[k+d \atop k\right] = \frac{c_{n,d,k}}{c_{n,0,k}} = \frac{\Gamma(\frac{k+d+1}{2})\Gamma(\frac{1}{2})}{\Gamma(\frac{k+1}{2})\Gamma(\frac{d+1}{2})}.
\]
Now, we have the formula for the expected Lipschitz--Killing curvature of the excursion set, which is proportional to the Minkowski functional.
\begin{theorem}
\label{thm:minkowski}
Under Assumptions \ref{as:B-R'} and \ref{as:A-T}, for $k=0,1,\ldots,n$,
\begin{align}
\label{minkowski}
 \mathbb{E}[\mathcal{L}_{k}(E_v)]
 =& \omega_{n-k}^{-1}\binom{n}{k}\mathbb{E}[\mathcal{M}_{n-k}(E_v)] \\
 =& \sum_{d=0}^{n-k} \left[k+d \atop k\right] \mathcal{L}_{k+d}(E)\,\Xi_{d}(v) + o(\nu^{s_0-2}). \nonumber
\end{align}
\end{theorem}

\begin{remark}
\label{rem:GKF}
The Gaussian Kinematic formula (GKF) established by \cite{taylor-adler:2009,adler-taylor:2007,adler-taylor:2019} states that the expected Euler characteristic for the Gaussian related field is written as a linear combination of the Lipschitz--Killing curvature $\mathcal{L}_d(E)$ multiplied by the Euler characteristic density $\rho_d(v)$ defined by the marginal distribution of the field.
Here, the Gaussian related field extended by \cite{panigrahi-etal:2019} is a random field defined by
\[
 X(t)=F(Z^{(1)}(t),\ldots,Z^{(p)}(t)), \quad t\in E,
\]
where $Z^{(1)},\ldots,Z^{(p)}$ are independent smooth random fields on $E$ with marginal mean 0 and variance 1.
The expression (\ref{minkowski}) is the same as the GKF formula, because $\gamma^{n/2}\mathcal{L}_d(E)$ is the Lipschitz--Killing curvature defined by the metric induced by the covariance $\mathrm{Cov}(\nabla X(t))=\gamma I_n$.
This similarity suggests the possibility of the extension of the GKF formula.
\end{remark}

\begin{remark}
\label{rem:contiguous}
The perturbation term $\Delta_{1,n}(x)$ in (\ref{D1n}) of Theorem \ref{thm:perturbation} determines the local power of the Minkowski functional for testing Gaussianity.
Suppose a setting that we have $N$ i.i.d.\ random fields $X^{(1)}(\cdot),\ldots,X^{(N)}(\cdot)$,
and let $E_v^{(i)}$ be the excursion set of $X^{(i)}(\cdot)$.
Define the test statistic
\[
 T_{N} = \frac{1}{\sqrt{N}}\sum_{i=1}^N\bigl(\mathcal{L}_k(E_v^{(i)})-\mu_{k,v}\bigr)/\sigma_{k,v},
\]
where $\mu_{k,v}=\mathbb{E}_0[\mathcal{L}_k(E_v)]$ and $\sigma_{k,v}^2=\mathrm{Var}_0(\mathcal{L}_k(E_v))$ are mean and variance under the Gaussian hypothesis $\nu=0$, respectively.
Let
\[
\delta_{k,v} = \sum_{d=0}^{n-k} \left[k+d \atop k\right] \mathcal{L}_{k+d}(E)\,
 \gamma^{-d/2} (2\pi)^{-d/2} \phi(v) \Delta_{1,d}(v).
\]
Then, under the contiguous alternative $\nu=c/\sqrt{N}$,
because $\mathrm{Var}(\mathcal{L}_k(E_v))=\mathrm{Var}_0(\mathcal{L}_k(E_v))(1+o(1))$,
the limiting distribution of $T_N$ when $N\to\infty$ is $N\bigl(c\delta_{k,v}/\sigma_{k,v},1)$.
\end{remark}

\subsection{A numerical example}
\label{sec:numerical}

As an example of a weakly non-Gaussian random field,
we consider a Gaussian random field by adding a quadratic form using an isotropic kernel.
By choosing a Gaussian kernel, the $N$-point correlation functions can be expressed explicitly.
We set the index set $E$ to be a square
$E=[0,e_1]\times\cdots\times[0,e_n]$
with the Lipschitz--Killing curvature
$\mathcal{L}_d = \prod_{i_1<\cdots<i_d} e_{i_1}\cdots e_{i_d}$.

First we define a random field $X(t)$ on $E$ with mean 0 and covariance function
\begin{equation}
\label{covg}
 \mathbb{E}[X(s) X(t)] = \rho(s,t) = e^{-\frac{g}{2}\Vert s-t\Vert^2}.
\end{equation}
We define
\begin{equation}
\label{delta-tau}
 Y(t) = X(t) + \delta \int_{\mathbb{R}^n} X(t) X(t') h(t',t) \dd t', \quad
 h(s,t) = \frac{1}{(2\pi\tau)^{n/2}}e^{-\frac{1}{2\tau}\Vert s-t\Vert^2},
\end{equation}
by the isotropic Gaussian kernel convolution, and let
\[
 Z(t) = (Y(t)-\mathbb{E}[Y(t)])/\omega_\delta, \quad \omega_\delta^2 = \mathrm{Var}(Y(x)).
\]

Noting that
\begin{align*}
 (\rho\ast h)(s,t)
=& \frac{1}{(1+g\tau)^{n/2}} e^{-\frac{1}{2}\frac{g}{1+g\tau}\Vert s-t\Vert^2}, \\
 (h\ast\rho\ast h)(s,t)
=& \frac{1}{(1+2g\tau)^{n/2}} e^{-\frac{1}{2}\frac{g}{1+2g\tau}\Vert s-t\Vert^2},
\end{align*}
we have
\begin{align*}
\mathbb{E}[Y(t)]
=& \delta \int \mathrm{cum}(X(t) X(t')) h(t',t) \dd t'
= \delta \int \rho(t,t') h(t',t) \dd t'
= \delta (\rho\ast h)(t,t),
\end{align*}
\begin{align*}
\mathrm{cum}(Y(s),Y(t))
=& \mathrm{cum}(X(s),X(t)) \\
&+ \delta^2 \int\int \mathrm{cum}(X(s) X(s'),X(t) X(t')) h(s',s) h(t',t) \dd s' \dd t'\\
=& \rho(s,t) + \delta^2 \int\int \rho(s,t)\rho(s',t') h(s',s) h(t',t) \dd s' \dd t'\\
&+ \delta^2 \int\int \rho(s,t')\rho(s',t) h(s',s) h(t',t) \dd s' \dd t'\\
=& \rho(s,t) + \delta^2  \rho(s,t) (h\ast\rho\ast h)(s,t)
 + \delta^2 (\rho\ast h)(s,t) (\rho\ast h)(t,s).
\end{align*}
The higher-order cumulants of $Y(t)$, and then those of $Z(t)$ can be obtained similarly.
Finally we have the 2, 3, 4-point correlation functions of $Z(t)$ as
\begin{align*}
\rho(x)
=& \biggl[e^{-g x}
 + \delta^2 \frac{1}{(1+2g\tau)^{n/2}} e^{-2g\frac{1+g\tau}{1+2g\tau} x}
 + \delta^2 \frac{1}{(1+g\tau)^{n}} e^{-\frac{2g}{1+g\tau} x}\biggr]\Big/\omega_\delta^2,
\end{align*}
where
\[
 \omega_\delta^2
 = 1
  + \delta^2 \frac{1}{(1+2g\tau)^{n/2}}
  + \delta^2 \frac{1}{(1+g\tau)^{n}},
\]
\begin{align*}
\kappa^{(3)}(x_1,x_2,x_3)
=& \biggl[\delta \sum_{i_1\ne i_2} \frac{1}{(1+g\tau)^{n/2}}
 e^{-g x_{i_1}}
 e^{-\frac{g}{1+g\tau} x_{i_2}} \\
 &+ \delta^3 \sum_{(i_1,i_2,i_3)\in S_3} \frac{1}{(1+g\tau)^{n/2}}\frac{1}{(1+2g\tau)^{n/2}}
 e^{-g x_{i_1} -\frac{g}{1+g\tau} x_{i_2} -\frac{g}{1+2g\tau} x_{i_3}} \\
 &+ 2\delta^3 \frac{1}{(1+g\tau)^{3n/2}}
 e^{-\frac{g}{1+g\tau}(x_1+x_2+x_3)}\biggr]\Big/\omega_\delta^3,
\end{align*}
where $S_3$ is the set of permutations of $\{1,2,3\}$ ($|S_3|=6$), and
\begin{align*}
\kappa^{(4)}(x_{(12)},x_{(13)},\ldots,x_{(34)})
= \biggl[
 2 \delta^2 \sum_{(e_1,e_2,e_3)\in DC_4} \frac{1}{(1+g\tau)^{n}}
 e^{-g x_{e_1}}
 e^{-\frac{g}{1+g\tau}(x_{e_2}+x_{e_3})} \\
+ 16 \delta^4 \sum_{(e_1,e_2,e_3,e_4)\in L_4} \frac{1}{(1+g\tau)^{2n}} e^{-\frac{g}{1+g\tau}(x_{e_1}+x_{e_2}+x_{e_3}+x_{e_4})}
\biggr]\Big/\omega_\delta^4,
\end{align*}
where
$DC_4 = \{ (e_1,e_2,e_3) :$ directed chain connecting 4 vertices $1,2,3,4$ of length 3$\}$ ($|DC_4|=24$), and
$L_4 = \{ (e_1,e_2,e_3,e_4) :$ undirected loop passing through 4 vertices $1,2,3,4$ of length 4$\}$ ($|L_4|=3$).
For example, $((13),(23),(24))\in DC_4$ and $((13),(23),(24),(14))\in L_4$.
From these formulas, the derivatives of the $N$-point correlation functions at the origin are evaluated explicitly.

\begin{figure}
\begin{tabular}{cc}
\scalebox{0.35}{\includegraphics{./ec-50-0-0-200-100-4.pdf}}
&
\scalebox{0.35}{\includegraphics{./ec-50-0_5-0_1-200-100-4.pdf}}
\end{tabular}
\caption{Euler characteristic curves for a Gaussian random field (left) and a non-Gaussian random field (right). \newline (solid: simulation, dotted: no correction, dot-dashed: skewness, dashed: skewness+kurtosis)}
\label{fig:convolution}
\end{figure}

Figure \ref{fig:convolution} depicts the simulated Euler characteristic curves (solid line) and the expected Euler characteristic curves without correction (dotted line), with corrections using the 3-point correlation $\kappa^{(3)}$ (dot-dashed line), and with corrections using 3- and 4-point correlations $\kappa^{(3)}$, $\kappa^{(4)}$ (dashed line).
In this study, we set the 2-dimensional index set as $E=[0,1]^2$, and the parameter in (\ref{covg}) $g=50$.
For the non-Gaussian case, we set as $\delta=0.5$ and $\tau=0.1$ in (\ref{delta-tau}).
The square index set $E=[0,1]^2$ is divided into $N\times N$ parts ($N=200$), and the data is generated
at the lattice points $\{(i_1/N,i_2/N) \mid i_1,i_2\in\mathbb{Z} \}$ by the R package \texttt{RandomFields}.
For the convolution, the R package \texttt{imagine} was used.

We approximate the excursion set by a triangle complex as in Figure \ref{fig:excursion} (right) and compute the approximate Euler characteristic.
The Euler characteristic for the threshold $v$ is
\[
 \chi(v) = \#\mbox{vertices}(v) - \#\mbox{edges}(v) + \#\mbox{triangles}(v).
\]
To calculate the Euler characteristic curve, we start form the initial value $\chi(v)=0$ when $v=\infty$, and then change the value of $v$ from $\infty$ to $-\infty$.
$\chi(v)$ changes at a critical threshold $v=v^*$ satisfying $X(t)=v^*$ at a vertex $t$.
Then, we calculate the difference $\chi(v^*+0)-\chi(v^*)$ by counting the numbers of the appearing/disappearing edges and triangles, all of which are connected to the newly appearing vertex at $v^*$.
They are easy to identify.

This procedure is based on the same principle as the standard proof of Morse's theory \cite{milnor:1963}.
The extension to 3 or higher dimensional case is straightforward.

\section{Cosmic data analysis: Future prospects}
\label{sec:cosmic}

Finally, we discuss briefly how the obtained formula will be utilized in the context of cosmic research.

In cosmology, the research is often based on data by simulator.
In order to check the validity of a cosmic model, the Minkowski functional curves are calculated based on simulation data generated by the assumed cosmic model, and are compared with the observed Minkowski functional curves to see the goodness-of-fit.
If the goodness-of-fit is not good, we conclude that the assumed model is not true.
Because the computational cost of the simulator is very high, the perturbation formulas will be very helpful in research if it provides the same information as simulators.

The Euler characteristic density $\Xi_n(v)$ in (\ref{Xin}) when the dimension is $n=2$ was derived and examined by \cite{matsubara:2010}.
To check the accuracy of the formula $\Xi_2(v)$,
he considered a simple model for the temperature map of the CMB radiation:
The weak non-Gaussian temperature fluctuations in the direction $\theta$ (spherical coordinates), denoted by $X(\theta)=(T(\theta)-\mathbb{E}[T(\theta)])/\mathbb{E}[T(\theta)]$, are generated by
\begin{equation*}
 X(\theta) = X_\mathrm{G}(\theta) - 3 f_\mathrm{NL} X_\mathrm{G}(\theta)^2 + 9 g_\mathrm{NL} X_\mathrm{G}(\theta)^3 - (\mathrm{monopole} + \mathrm{dipole}),
\end{equation*}
where $X_\mathrm{G}(\theta)=(T_\mathrm{G}(\theta)-\mathbb{E}[T_\mathrm{G}(\theta)])/\mathbb{E}[T_\mathrm{G}(\theta)]$ is the auxiliary temperature fluctuation of a Gaussian random field with a Saches--Wolfe power spectrum.
The local relation between $X_\mathrm{G}$ and $X$ is a simple assumption of the toy model just for the numerical demonstration, and is not necessarily accurate in reality.
This toy model is derived from the following assumptions: the primordial non-Gaussianity is given by a local-type ansatz with parameters $f_\mathrm{NL}$ and $g_\mathrm{NL}$, and the CMB temperature fluctuations are given by an asymptotic formula of the Sachs--Wolfe effects on large-angular scales beyond the horizon at recombination epoch.
Then, a radial Gaussian filter is applied to $X(\theta)$.

In this setting, 100,000 realizations of the non-Gaussian CMB temperature map are generated with random seeds for the Gaussian random field $X_\mathrm{G}(\theta)$, and Minkowski functionals for filtered $X(\theta)$ are calculated in each realizations, and finally they are averaged over.
The numerical results are compared with the theoretical curves as functions of the threshold.
It is shown that the formula for Minkowski functionals $\mathcal{M}_0$, $\mathcal{M}_1$ and $\mathcal{M}_2$ are extremely accurate within small sample variances among the 100,000 realizations (Figure 1 of \cite{matsubara:2010}).
In general, the CMB temperature fluctuations with a primordial non-Gaussianity is not given by the simply local relation of the toy model discussed above.
They can have much complicated higher-order cumulants.
Our formula is applicable to these general cases of weakly non-Gaussian fields with arbitrary forms of cumulants.

In this paper, we have obtained the Minkowski functionals in arbitrary dimensions $n$, and on a bounded domain.
In cosmology, the cases $n=1,2,3$ are interesting, because we have many cosmological fields that are obtained by various kinds of survey observations, such as 2-dimensional maps of the CMB, weak-lensing field of galaxy surveys, and 3-dimensional maps of the large-scale structure of the universe probed by galaxy distributions.
We are preparing a separate paper on the 3-dimensional Minkowski functionals, which compares the Minkowski functionals in cosmological simulations of a large-scale structure with theoretical predictions of the weakly non-Gaussian formula, combined with the theoretical predictions of non-Gaussianity parameters by nonlinear perturbation theory \cite{matsubara-hikage-kuriki:2020}.

The boundary effects of the Minkowski functionals are also of interest in cosmological applications, because there are always survey boundaries in the actual data.
The formula of Theorem \ref{thm:boundary}
may find applications in cosmological fields with analytic corrections of the boundary effects.
The advantage of the boundary corrections may be significant in repetitive small angle surveys to reduce the effect of the so-called cosmic variance.


\section{Introduction}
\vspace{-.1in}
Cable-driven robots are gaining increasing attention due to their adaptiveness and safety. Tensegrity structures have many applications: from manipulation \cite{lessard2016bio}, locomotion \cite{sabelhaus2018design}, morphing airfoil \cite{chen2020design} to spacecraft lander \cite{bruce2014superball}. While useful and versatile, they are difficult to accurately model and control. Identifying system parameters is necessary, either to learn controllers in simulation (as real-world experiments are time-consuming, expensive and dangerous), or for traditional model-based control. In all these cases, the spring-rod representation considered in this work is the basic modeling element.

However, the spring-rod system has high degrees of freedom for system identification. Physics-based methods for simulation require accurate models that capture non-linear material behavior, which are difficult to construct. In contrast, data-driven methods can simulate any system from observed data, with sufficient training data. But the large number of variables and non-linear material properties necessitate copious amounts of training data.

Motivated by these issues, we propose a data-driven differentiable physics engine that combines the benefits of data-driven and physics-based models, while alleviating most of their drawbacks, and is designed from first principles. Previous data-driven models have required large amounts of data, because they learn the  parameters \emph{and} the physics of the system.  Furthermore, the hidden variables and black box nature of these models are not explainable, and difficult to transfer to new environments. Our approach is based on the observation that the equations that govern the motion of such systems are well-understood, and can be directly baked into the data-driven model. Such a design can reduce demands on training data and can also generalize to new environments, as the governing principles remain the same. We further simplify the differentiable engine by modularization, which compartmentalizes the problem of learning the dynamics of the whole system to smaller well-contained problems. For each module, we also reduce the dimension from 3D to 1D, by taking advantage the properties of spring-rod systems, which allows for efficient parameter inference using linear regression. As a side benefit, the regression parameters correspond to physical quantities, such as the spring stiffness or the mass of the rod, making the framework explainable.

\begin{figure}
\vspace{-0.7cm}
    \centering
    \includegraphics[height=0.25\linewidth, trim={6.5cm 5cm 2cm 1cm}, clip]{figures/simple_spring_rod_system.pdf}
    \includegraphics[height=0.25\linewidth, trim={0cm 0cm 0cm 0cm}, clip]{figures/super_ball_ground.jpg}
    \includegraphics[height=0.25\linewidth]{figures/tensegrity_real_platform.jpg}
    \vspace{-2mm}
    \caption{A basic element with one rod connected by two springs (left). A complex assembly of rods and springs forming a tensegrity robot in simulation (middle), and the real world (right).}
    \label{fig:spring_rod_model}

\end{figure}

\vspace{-3mm}
\section{Related Work}
\vspace{-.1in}
Traditional methods for system identification build a dynamics model by minimizing the prediction error~\cite{swevers1997optimal} \cite{hansen2001completely}.
These methods require parameter refinement and data sampling in an iterative fashion, to decrease the prediction error. This iterative process can be avoided using data-driven techniques that directly fit a physics model to data~\cite{rosenblatt1958perceptron,rumelhart1986learning,asenov2019vid2param}. However, these techniques treat the dynamics as a black box, are data hungry, and require retraining in a new environment.

Instead of treating the environment as a black box, \emph{interaction network} \cite{battaglia2016interaction} took the first step to modularize objects and their interactions. Later, researchers extended this idea by a hierarchical relation network for graph-based object representation~\cite{mrowca2018flexible} and a propagation network \cite{li2019propagation} for the multi-step propagation network. While these methods are an improvement over previous approaches, they still treat the interactions between different objects as black boxes and try to learn them from data, even though the governing equations of motion are well-understood.

Quite a few authors have recently introduced differentiable physics engines that focus on many aspects not central to our work. For example, forward dynamics ~\cite{heiden2019interactive}, Material Point Method (MPM) \cite{hu2019chainqueen}, linear complementarity problems (LCP)~\cite{de2018end} , augmented Lagrangian~\cite{landry2019differentiable}, differentiable programming~\cite{hu2019difftaichi}, augmented traditional physics simulators\cite{sain}, and LSTM dynamics model without system topology~\cite{golemo2018sim}. Researchers have also proposed differentiable engines specific to certain kinds of objects, such as molecules~\cite{jaxmd2019}, fluids~\cite{spnets2018}, and cloth~\cite{liang2019differentiable}. Recent works on tensegrity robots~\cite{surovik2019adaptive, littlefield2019kinodynamic, surovid2018any} make major improvement on locomotion in simulation and have great challenges on policy transfer to real world system. All these works motivate us to mitigate the reality gap for cable driven robots between simulation and reality via system identification.

\vspace{-3mm}
\section{Methods}
\vspace{-0.1in}

\begin{figure}
\vspace{-0.3cm}
    \centering
    \includegraphics[width=\textwidth]{figures/simple_physics_nn_flow_chart.pdf}
    \vspace{-7mm}
    \caption{Flow chart showing the data flow when simulating one time step with our physics engine.}
    \label{fig:architecture}
    \vspace{-5mm}
\end{figure}

Our system views a spring-rod system as a composition of basic \emph{elements} (see Fig.~\ref{fig:spring_rod_model}(left)), where springs generate forces that influence rod dynamics. We subdivide each time step of the simulation into three modules: force generation, acceleration computation, and state update/integration (see Fig.~\ref{fig:architecture}). The physics engine takes as input the current rod state $S_t = \{\boldsymbol{p, v, q, \omega}\}$, where $p$ is position, $v$ is linear velocity, $q$ is orientation (expressed as a quaternion), and $\omega$ is the angular velocity. Based on $S_t$, the position and linear velocity $p^{m}, v^{m}$ of the two rod endpoints is computed, and is used to compute the relative compression (or expansion) and velocity $\Delta p^{m}, \Delta v^{m}$ of the two attached springs. Then, the first module predicts the spring forces $f$, the second module computes the linear and angular accelerations $a,\alpha$, and the third module integrates the new state $S_{t+1}$.

\vspace{-2mm}
\subsection{System Topology Graph}
\vspace{-2mm}

\begin{wrapfigure}[4]{r}{0.38\textwidth}
\vspace*{-6mm}
\begin{overpic}[width=0.38\textwidth]{figures/topology_graph.pdf}
    \put(1.5,8){{\small spring1}}
    \put(26,13){{$f_1$}}
    \put(45,8){{\small Rod}}
    \put(69.5,13){{$f_2$}}
    \put(80.5,8){{\small spring2}}
\end{overpic}
\vspace{-5mm}
\caption{\small Element topology graph.}
\label{fig:topology}
\end{wrapfigure}

We use a topology graph to represent interconnections between different components of the spring-rod system. Each rod and spring has a corresponding vertex, and directed edges represent relations between them. Figure~\ref{fig:topology} shows an example topology graph for the basic spring-rod element shown in Figure~\ref{fig:spring_rod_model}(left).

\vspace{-3mm}
\subsection{Force Generation Module}
\vspace{-3mm}

The relative compression (or expansion) $\Delta p^{m}$ and velocity $\Delta v^{m}$ of each spring is given as input to the force generation module, which outputs the spring forces $f$  by Hooke's law. As shown in Fig.\ref{fig:spring_physics_engine}, two unknown parameters, stiffness $K$ and damping $c$, can be easily learned using linear regression.

\begin{figure
    \centering
    \includegraphics[width=0.8\textwidth]{figures/spring_physics_engine_graph.jpg}
    \vspace{-1mm}
    \caption{Force generation module, which uses dimensionality reduction to compute spring forces.}
    \label{fig:spring_physics_engine}
 \vspace{-5mm}
\end{figure}

\vspace{-3mm}
\subsection{Acceleration Generation Module}
\vspace{-3mm}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{figures/rod_physics_engine_graph_updated.jpg}
    \vspace{-2mm}
    \caption{Acceleration generation module}
    \label{fig:rod_physics_engine}
    \vspace{-6mm}
\end{figure}
The acceleration generation module takes the spring forces $f$ and control force $f_u$ as input and outputs each rod's linear and angular accelerations $a,\alpha$ as shown in Fig~\ref{fig:rod_physics_engine}. $f_1$ and $f_2$ are spring forces on the two rod ends, $f_u$ is control force, $r$ is the half-length rod vector, $r_u$ is control force arm, $R$ is the rod local/world frame rotation matrix. The rod mass $M$ and inertia $I$ are unknown parameters to identify.

\vspace{-3mm}
\subsection{Integration Module and Method Implementation}
\label{integraion_module}
\vspace{-3mm}

The integration module computes forward dynamics of each rod using the current accelerations $a,\alpha$. We apply the semi-implicit Euler method~\cite{stewart2000implicit} to compute the updated state $S_{t+1}=\{p_{t+1}, v_{t+1}, q_{t+1}, \omega_{t+1}\}$ at the end of the current time step.

The learning module receives the current state $S_t$ and returns a prediction $\hat{S_{t+1}}$.  The loss function is the MSE between the predicted $\hat{S_{t+1}}$ and ground truth state $S_{t+1}$. The proposed decomposition, first-principles approach and the cable’s linear nature allow the application of linear regression, which helps with data efficiency. This linear regression step has been implemented as a single layer neural network without activation function on pyTorch \cite{paszke2019pytorch}.

\vspace{-3mm}
\section{Experiments}
\vspace{-3mm}
The task is to identify parameters including spring stiffness $K$, damping $c$ and rod mass $M$, inertia $I$.
\vspace{-3mm}
\subsection{Simple Spring-Rod Model Identification}

\begin{figure}
\vspace{-0.3cm}
    \centering
    \begin{minipage}[b]{.32\linewidth}
    \includegraphics[height=1.4cm, trim={1.8cm 0cm 0cm 0cm}, clip]{figures/interaction_nn_arch.png}
    \vspace{-6mm}
    \caption{Interaction Network}
    \label{fig:interaction_network_architecture}
    \end{minipage}%
    \hfill
    \begin{minipage}[b]{.67\linewidth}

    \includegraphics[width=\textwidth]{figures/koopman_w_integration.jpg}
    \vspace{-6mm}
    \caption{Koopman with Integration Module}
    \label{fig:koopman_w_int_architecture}
    \label{fig:interaction_network_w_int_architecture}
    \end{minipage}

\vspace{-0.2in}
\end{figure}
First, we identify these parameters in a simple spring-rod system as shown in Figure~\ref{fig:spring_rod_model} (a). \textbf{Interaction} is an improved version of the Interaction Network \cite{battaglia2016interaction} as shown in Fig. \ref{fig:interaction_network_architecture}. It has two Multilayer Perceptrons (MLPs), one to generate spring forces $f$ and the other to generate rod state $S_{t+1}$. Unlike \cite{battaglia2016interaction}, which takes raw state $S_t$ as input, we generate $\hat{\Delta p^{m}_t}, \hat{\Delta v^{m}_t}$ as input. \textbf{Interaction+Int} appends the integration module to the Interaction Network, and replaces input $S_t$ of $MLP_2$ by $r$.  \textbf{Koopman+Int} is to use the Koopman operator to predict accelerations and apply the Integration Module to map them to $S_{t+1}$.  \textbf{Interaction} only predicts a $S_{t+1}$ in training data that is close to $S_t$. \textbf{Interaction+Int} experiences increasing error from accumulated prediction errors. The Koopman operator \textbf{Koopman+Int} designed from first principles gives accurate predictions similar to \textbf{Ours} in this simple system. Comparison of errors is shown in Figure~\ref{fig:superball_system_identification} a).

\begin{figure}

    \begin{minipage}[t]{0.5\linewidth}

    \centering
    \includegraphics[height=0.41\linewidth]{figures/single_rod_known_topology_position_error.pdf}
    \includegraphics[height=0.41\linewidth]{figures/single_rod_known_topology_quaternion_error.pdf}

    \caption*{a) Simple Spring-Rod Model Identification}
    \end{minipage}
    \begin{minipage}[t]{0.5\linewidth}

    \centering
    \includegraphics[height=0.41\linewidth]{figures/super_ball_same_params_position_error.pdf}
    \includegraphics[height=0.41\linewidth]{figures/super_ball_same_params_quaternion_error.pdf}

    \caption*{b) Complex Tensegrity Model Identification}
    \end{minipage}

    \caption{Physical Parameters Estimation}
    \label{fig:superball_system_identification}
    \vspace{-0.6cm}
\end{figure}

\vspace{-.1in}
\subsection{Complex Tensegrity Model Identification}
\vspace{-.1in}
We consider an icosahedron tensegrity system as shown in Fig. \ref{fig:spring_rod_model} (c). It is composed of 6 rods and 24 springs. Each rod is connected to 8 springs and has a length of 1.04m. Each spring's rest length is 0.637m. We set the gravity constant to $g=-9.81$ in Mujoco. We collect 1000 trajectories with different initial conditions for training, 200 for validation and 100 for testing. The result is shown in Figure \ref{fig:superball_system_identification} b). \textbf{Our} approach outperforms \textbf{Koopman+Int} because designing basis functions for The Koopman operator has an increased data requirement relative to our approach.

\vspace{-.1in}
\subsection{Data Efficiency Experiment}
\vspace{-.1in}
The proposed method has relatively small data requirements as shown in Fig. \ref{fig:superball_diff_params_system_data_efficient_and_generalization} a). Instead of training on 1000 trajectories, which have 736,167 time steps in total, we train our model with less data and evaluate performance. We randomly select 10\%, 1\%, 0.1\%, 0.01\%  of the 736,167 time steps for training. The model achieves good performance even with 73 time steps for training. All trajectories are from the complex tensegrity setup with different parameters.

\begin{figure}[ht]
\vspace{-.1in}
\begin{minipage}[t]{.5\linewidth}
\vspace{0pt}
\centering
\includegraphics[width=0.49\linewidth]{figures/super_ball_data_efficient_position_error.pdf}
\includegraphics[width=0.49\linewidth]{figures/super_ball_data_efficient_quaternion_error.pdf}
\vspace{-0.5cm}
\caption*{a) Efficiency on Various Training Set}
\end{minipage}%
\begin{minipage}[t]{.5\linewidth}
\vspace{0pt}
\centering
\includegraphics[width=0.49\linewidth]{figures/super_ball_w_control_force_generalization_long_traj_position_error.pdf}
\includegraphics[width=0.49\linewidth]{figures/super_ball_w_control_force_generalization_long_traj_quaternion_error.pdf}
\vspace{-0.5cm}
\caption*{b) Generalization on New Testing Set}
\end{minipage}
\vspace{-0.2cm}
\caption{Data Efficiency and Model Generalization Experiment.}
\label{fig:superball_diff_params_system_data_efficient_and_generalization}
\vspace{-0.15in}
\end{figure}

\textcolor{black}{The proposed solution achieves very low error at a magnitude of $10^{-5}$, since it: 1) introduces a first-principles approach in learning physical parameters; 2) removes redundant data from regression; 3) operates -for now- over relatively clean data from simulation before moving to real-world data. }

\vspace{-.1in}
\subsection{Model Generalization Experiment}
\vspace{-.1in}
\begin{wrapfigure}{r}{0.5\textwidth}
    \vspace{-.3in}
    \centering
    \includegraphics[width=0.49\textwidth]{figures/rod_control_force_scalar.jpg}
    \vspace{-.1in}
    \caption{Control Force Scalar}
   \label{fig:control_force_scalar}
 \vspace{-0.2in}
\end{wrapfigure}
This section generalizes the physics engine trained with a dataset without external forces to a dataset with such forces. We are interested in evaluating: 1) how the physics engine performs for longer time horizons (e.g., after 2000 time steps); 2) if it can adapt to new scenarios. We generate a new dataset with 20,000 time steps trajectories with random directed perturbation forces $f_u$. The external force $f_u$ does not have the same scale as the internal spring forces, so we add a new scalar module with only one parameter $h$, as in Fig. \ref{fig:control_force_scalar}.  We also apply dimensionality reduction to improve data efficiency. The tuning process is to freeze all other modules' weights and train this by the new dataset. The error graphs are shown in Fig. \ref{fig:superball_diff_params_system_data_efficient_and_generalization} b).
\vspace{-.1in}
\section{Conclusion and Future Work}
\vspace{-.1in}
This paper proposes a differentiable physics engine for system identification of spring-rod systems based on first principles. The engine has three modules: force generation, acceleration generation and integration, which express the corresponding physical processes of spring-rod systems. This results in reduced data requirements and improved parameter accuracy. It also provides an explainable, accurate and fast physics engine. In the future, we plan to address contacts and friction. This will involve replacing the linear regressor with nonlinear models in the existing modules. To overcome noise in real data, we plan the addition of a residual network along with the nonlinear model. These changes may also help with temporal scalability. The ultimate mission, we generate policy from our identified engine and evaluate on the real platform to finally mitigate the reality gap.

\small
\bibliographystyle{plain}