\documentclass[epub,pdf]{ebook}
\usepackage[cm]{sfmath}
\usepackage{amsmath}
\usepackage{color}
\usepackage[pdftoolbar=false, pdfpagemode=UseNone, linkbordercolor={0 0 0}]{hyperref}

\ident{ES}{19}{6}{2019}{2019ES000651}
\cpright{}{2019}
\pagestyle{empty}
\definecolor{brown}{RGB}{157,23,4}

\begin{document}

\vspace*{2\baselineskip}

\title{Kaula's rule -- a consequence\\ of probability laws by\\ A. N. Kolmogorov\\ and his school}

\author{E. B. Gledzer$^{1}$, G. S. Golitsyn$^{1}$}

\affil{}{A. M. Obukhov Institute of Atmospheric Physics RAS, Moscow,  Russia.}	

\abstract{The paper analyzes a possible cause for the universal behavior of the covariating fluctuations of planet gravity. The consideration based on the idea that the topography fluctuations are governed by a random Markov process leads to universal dependence $k^{-2} $with $k$ being the amplitudes of the $j$-th spherical harmonic of the terrain profile. This law known as the Kaula's rule is then derived from the solution of the Fokker -- Plank equation for the fluctuations of the terrain profile as the function of the horizontal coordinate. The respective diffusivities for Earth and Venus are retrieved from the existing \linebreak\newpage}\noindent\centerline{\parbox{100mm}{  data.}}


\section{Introduction}\vspace*{-.5pc}

The Kaula's rule [\itc{Kaula,} \reflink{kaul66}{1966}] sometimes called ``the rule of thumb'' exists for over half a century. However, its origin is unknown. In its original form it claims that the terrestrial gravity field fluctuations expanded in spherical harmonics have the spectral asymptotics proportional to $j^{-2}$, where $j$ is the number of the harmonic. Similar behavior was established to hold on Moon, Mars [\itc{Rexer and Hirt,} \reflink{rexe15}{2015}], and Venus [\itc{Turcotte,} \reflink{turc97}{1997}]. In this decade the same property has been found for asteroid Vesta whose diameter does not exceed 300 km [\itc{Konopliv et al.,} \reflink{kono14}{2014}] and for asteroids of few kilometers in size [\itc{McMahon et al.,} \reflink{mcma16}{2016}]. In later 1980-s it was found that surface reliefs of Earth and Venus have the same $j^{-2}$ asymptotics [\itc{Turcotte,} \reflink{turc97}{1997}]. In [\itc{Rexer and Hirt}, 2015] the harmonics with $j$ up to tens thousands have been computed. This means that gravity field fluctuations and relief are connected. Therefore ``the rule of thumb'' is expected to be a law whose nature still awaiting for explanation.

Here we use the Fokker -- Planck equation in the form presented by \itc{A.~N. Kolmogorov,} [\reflink{kolm34}{1934}] and its further results of its statistical analysis developed by his pupils A.~M. Obukhov [\itc{Obukhov,} \reflink{obuk59}{1959}] and by \itc{A.~S. Monin and A.~M. Yaglom,} \reflink{moni75}{1975}]. The basic equation formulated in [\itc{Kolmogorov,} \reflink{kolm34}{1934}] is more general than it is necessary here. Later this equation was used for studying other natural phenomena: like turbulence in [\itc{Obukhov,} \reflink{obuk59}{1959}], and recently sea wind waves spectra, the cumulative distribution of litospheric plates, etc. in [\itc{Golitsyn,} \reflink{goli18}{2018}]. This impelled us to return to an attempt to derive analytically the Kaula's rule. For reader's convenience we present the description of ideas and methods of  [\itc{Kolmogorov,} \reflink{kolm34}{1934}; \itc{Obukhov,} \reflink{obuk59}{1959}; \itc{Monin and Yaglom,} \reflink{moni75}{1975}] rarely used in general geophysics but still appearing in various papers and books during past decades.\vspace*{-.5pc}


\section{Basic Equations}\vspace*{-.5pc}


In the beginning of 1930-s A.~N. Kolmogorov had published three papers on analytical methods in the probability theory. The two-page work [\itc{Kolmogorov,} \reflink{kolm34}{1934}] contains the essence of his approach started by A. Einstein and developed further by Fokker and Planck. Kolmogorov proposed a fundamental solution to the evolution equation for the probability distribution $p(u_{i}, x_{i}, t)$ of the 6D vector $u, x$ provided a delta correlated force causes the respective Markov process. This is an approximation of the processes where correlation time of the random forces is much shorter than the reaction time of the considered system. In 1D approximation the action of respective random accelerations $a_t$ with $\left\langle a(t_{1} )a(t_{2} )\right\rangle = a^{2} \delta (t_{1} - t_{2})$ can be described by the simplest equations [\itc{Gledzer and Golitsyn,} \reflink{gled10}{2010}],

\begin{equation} \label{1}
\frac{du(t)}{dt} = a(t), \frac{dx(t)}{dt} = u(t),
\end{equation}

\noindent where $u(t)$ and $x(t)$ are the velocity and displacement of a particle under random acceleration $a(t)$.

A.~N. Kolmogorov used a Fokker-Planck type equation [\itc{Kolmogorov,} \reflink{kolm34}{1934}] to describe the evolution of the probability density $p(u_{i}, x_{i}, t)$,

\begin{equation} \label{2}
\frac{\partial p}{\partial t} + u_{i} \frac{\partial p}{\partial x_{i}} =
D\frac{\partial ^{2} p}{\partial u_{i}^{2}} ,
\end{equation}

\noindent where $D$ is a coefficient of diffusion in the velocity space. The fundamental solution\footnote{ In the original paper [$  $Kolmogorov$  $, 1934] the middle term under exponent was absent.} to eq. \eqref{2} has the form (see also [\itc{Monin and Yaglom,} \reflink{moni75}{1975}], {\S}24.4):

\begin{equation} \label{3}
p(u_{i} ,x_{i} ,t) = \left(\frac{\sqrt{3}}{2\pi Dt^{2}} \right)^{3}
\exp \left[-\left(\frac{u_{i}^{2}}{Dt} - \frac{3u_{i} x_{i}}{Dt^{2}} + \frac{3x_{i}^{2}}{Dt^{3}} \right)\right].
\end{equation}

A.~M. Obukhov [\itc{Obukhov,} \reflink{obuk59}{1959}] was the first who analyzed and used this equation. He showed that the coefficient $D$ is proportional to the generation/dissipation rate of the kinetic energy of the random motions $D = \varepsilon /2$. The details of his analysis are presented in [\itc{Monin and Yaglom,} \reflink{moni75}{1975}; \itc{Golitsyn,} \reflink{goli18}{2018}]. The solution \eqref{3} demonstrates that the seeked probability distribution has the normal character. This solution has two scales (angular brackets mean ensemble averaging):

\begin{equation} \label{4}
\left\langle u_{i}^{2} \right\rangle = \varepsilon t, \left\langle x^{2} \right\rangle \equiv r = \varepsilon t^{3}, \varepsilon = 2D.  %%      \eqref{GrindEQ__4_}
\end{equation}

Combining the scales for velocity and coordinates Obukhov formulated the scale for the mixed moment $\left\langle u_{i} x_{i} \right\rangle =K$ in terms of $\varepsilon$ as

%%\begin{equation}\label{4'}
\centerline{\hspace*{10pc} $K=\varepsilon t^{2} .$ \hspace*{6pc} (4') \hspace*{1pc}}
%%\end{equation}

\begin{figure}                      % Fig 1
\figurewidth{90mm}
\setimage{}{}{115mm}{}{2019es000651-f01}
\fcaplabel{Temporal behaviour of the second velocity and coordinate moments for different particle ensemble sizes as in [\itc{Gledzer and Golitsyn,} \reflink{gled10}{2010}] (see text).}\vspace*{1pc}
\end{figure}


The time dependencies eqs \eqref{4} and (4') were verified numerically for the ensembles of $N$ randomly moving particles [\itc{Gledzer and Golitsyn}, \reflink{gled10}{2010}]. \figref{1} demonstrates that even for $N$ = 10 the relationships \eqref{4} fulfil satisfactorily as it is shown in \figref{1}. Moreover, on introducing the dimensionless variables to non-dimensional ones (denoted with tilde overhead)
$$t=\tilde{t}\tau , u=\tilde{u}(\varepsilon \tau )^{1/2} , r=\tilde{r}(\varepsilon \tau ^{3} )^{1/2} $$
reduces eq. \eqref{2} to the universal or selfsimilar form [\itc{Gledzer and Golitsyn,} \reflink{gled10}{2010}]

\begin{equation} \label{5}
\frac{\partial p}{\partial \tilde{t}} +\tilde{u}\frac{\partial p}{\partial \tilde{r}} =\frac{1}{2} \frac{\partial ^{2} p}{\partial \tilde{u}^{2} } .
\end{equation}

The scales \eqref{4} can be interpreted as the structure functions with zero initial conditions which allow us to introduce their spectral forms [\itc{Monin and Yaglom,} \reflink{moni75}{1975}]. If we replace $t$ by $r$ (see eq. \eqref{4}), exclude time from the second scale \eqref{4} and put it into the first scales and to (4') and substitute the result into eq. (4'), we obtain the results of Kolmogorov -- Obukhov of 1941 $D_{u}(r) = (\varepsilon r)^{2/3} $ and Richardson -- Obukhov [\itc{Obukhov,} \reflink{obuk41}{1941}; \reflink{obuk59}{1959}; \itc{Richardson,} \reflink{rich26}{1926}; \reflink{rich29}{1929}] for the eddy diffusion coefficient:

\vspace*{-.5pc}
\begin{equation} \label{6}
K(r)=\varepsilon ^{1/3} r^{4/3} .
\end{equation}

These are the results for the inertial interval of turbulence. The results as relations \eqref{4} are valid irrespective of the space dimension and the numerics prove it [\itc{Gledzer and Golitsyn,} \reflink{gled10}{2010}]. It is surprising that since 1920s nobody tried to understand why the Richardson's ``law of -- 4/3'' holds up to 2--3 thousands km in the horizontal direction [\itc{Richardson,} \reflink{rich26}{1926};  \reflink{rich29}{1929}], i.~e. the vertical component and full velocity isotropy are not important. In 1999 Lindborg [\itc{Lindborg}, \reflink{lind08}{2008}] analyzed the experimental data of commercial flights in terms of structure functions. He showed that the dependence $r^{2/3} $ is traced up to distances of 2 thousands km. His results in terms of the eddy diffusion explain the Richardson results of 1920-s. These facts can be explained starting with a two-page note of Kolmogorov, 1934 [\itc{Kolmogorov,} \reflink{kolm34}{1934}], and the following development of his results by \itc{Obukhov,} [\reflink{obuk59}{1959}] and \itc{Monin and Yaglom,} [\reflink{moni75}{1975}]. Only the Markov properties of the process are important rather than the full isotropy except in horizontal direction. Our numerics of 2010 [\itc{Gledzer and Golitsyn,} \reflink{gled10}{2010}] confirmed this independence of the space dimensionality for the second moments with the accuracy of the numerical coefficients up to 0 \eqref{1}.

\section{Spectral Representations}

For further development we need an analytical connection between the second moments and their spectral representations when both of them are power laws. Let us consider first the structure function $D(t) = At^{\gamma }$, $0 < \gamma < 2$. The relation to the spectrum is [\itc{Monin and Yaglom,} \reflink{moni75}{1975}]

\begin{equation} \label{7}
D(t)=2\int _{0}^{\infty }(1-\cos \omega t)E(\omega )d\omega  ,
\end{equation}

The inversion gives for the spectrum [\itc{Monin and Yaglom,} \reflink{moni75}{1975}; \itc{Yaglom,} \reflink{yagl55}{1955}; \reflink{yagl86}{1986}]:

\begin{equation} \label{8}
E(\omega ) = \frac{C}{\omega ^{\gamma +1} } , C=\frac{A}{\pi } \Gamma (\gamma +1)\sin \frac{\pi \gamma }{2} ,
\end{equation}

\noindent where $\Gamma$ is the gamma function. The length scale related to the mean square displacement of a randomly moving particle has therefore the spectrum $C\omega ^{-4}, C = A/\pi$  and for the velocity scale $\left\langle u^{2} \right\rangle = \varepsilon t$ the spectrum is $\varepsilon \omega^{-2}$. According to the terminology introduced by \itc{Yaglom,} [\reflink{yagl55}{1955}; \reflink{yagl86}{1986}] the spectrum $\omega ^{-2}$ corresponds to a process with stationary increments of the first order ($0 < \gamma < 2$), while the spectrums $\omega^{-4}$ has the stationary increments of the second order ($2 < \gamma < 4$). In order to avoid the singularity at $\omega \to 0$ the transformation kernel in eq. \eqref{7} should be $(1 -cos \omega t)^n$, with $n = 1$ for velocities, and $n = 2$ for coordinates.



\section{Toward Thumb Rule}

When a flying altimeter measures the relief, we obtain a temporal variations of height signal $h(t)$. The time $t$ here is related to horizontal coordinate as $y = ut$, $u$ being the flight velocity. The standard Fokker -- Planck equation with stochastic velocity fluctuations in this case is

$$\frac{\partial p}{\partial t} = D\frac{\partial ^{2} p}{\partial h^{2} } $$
If $y = ut$ it transforms to

\begin{equation} \label{9}
\frac{\partial p}{\partial y} =D_{1} \frac{\partial ^{2} p}{\partial h^{2} } , D_{1} =\frac{D}{u} .
\end{equation}

The diffusion coefficient $D_{1}$ has the dimension of length. But if the vertical $L_{z}$ and horizontal $L_{y}$ scales are different, then $D_{1} = L_{z}^{2} /L_{y}$ which reflexes clearly the diffusion character of the relief formation process during millions years. Equation \eqref{9} yields the scale and the structure functions for zero initial condition in the form:
$$\left\langle h^{2} (y)\right\rangle = D_{1} y$$
with the corresponding spectrum

\begin{equation}\label{10}
S(k) = D_{1}k^{-2}, k = 2\pi / \lambda_{y},
\end{equation}

\noindent $\lambda_{y}$ being the corresponding horizontal wave length. This equations applicable for small scale continuous areas. The spectrum of the relief slopes in this approximation is $S_{\zeta }(k) = k^{2}S(k) = D = const.$ The constancy of a spectrum is referred to as ``white noise''. The slope of the relief

\begin{equation} \label{11}
\frac{dh(y)}{dy} =\varsigma (y)
\end{equation}

\noindent for which we accept the hypothesis of $\delta$-correlated horizontal coordinates
\begin{equation} \label{12}
\left\langle \varsigma (y_{1} )\varsigma (y_{2} )\right\rangle = \varsigma ^{2} \delta (y_{1} -y_{2} )
\end{equation}

Equation \eqref{12} corresponds to the Markov character of forces in [\itc{Kaula,} \reflink{kaul66}{1966}; \itc{Rexer and Hirt,} \reflink{rexe15}{2015}; \itc{Turcotte,} \reflink{turc97}{1997}]. Here $\varsigma ^{2}$ may be interpreted as the mean square slope angle, a dimensionless value.

Now we return to the reliefs of concrete planets. Recently an extensive paper on ultrahigh spherical harmonics analysis for topography of Earth, Mars and Moon has appeared [\itc{Rexer and Hirt,} \reflink{rexe15}{2015}] with the horizontal resolution up to hundreds of meters. Much earlier analyses of topography for Earth and Venus were given in [\itc{Turcotte,} \reflink{turc97}{1997}]. In that book there are also the data for mountain, hilly and plain areas of hundred kilometers in size at the Oregon state. The flights were done in various directions in each type of terrain. For records in the range 1 -- 60 km the mean exponent value of these spectra were found [\itc{Golitsyn,} \reflink{goli03}{2003}] to be 2.03 $\pm$ 0.04. Unfortunately, there were no data in [\itc{Turcotte,} \reflink{turc97}{1997}] for spectral amplitudes. Obviously $k^{-2}$ small size spectrum is a continuation of what Kaula found for large scale harmonics.

In the book [\itc{Turcotte,} \reflink{turc97}{1997}] the spectral exponent 2 (and other statistical measures, such as Hausdorff measure, the fractal index $d$) were retrieved from the original records. But no explanation of particular values was proposed.

The analysis done by Kolmogorov [\itc{Kolmogorov,} \reflink{kolm34}{1934}; \itc{Monin and Yaglom,} \reflink{moni75}{1975}] several decades earlier was clearer and simpler. The only assumption on $\delta$-correlation of slopes applies for the processes in consideration, i.e. the Markov's hypothesis. Such assumption is used in many branches of theoretical physics meaning only that the system reaction time is much larger than the correlation time of stochastic influence on it [\itc{Golitsyn,} \reflink{goli03}{2003}].

\section{Spherical Harmonics}

For large planetary scale features the spherical harmonic analysis is used. Let us consider the relief on a sphere $z(\varphi, \theta)$ with a radius $r$, where $\varphi$ is longitude, $0\le \varphi < 2\pi$, $\theta$ is a co-latitude, $0\le \theta \le \pi$. Then
$$\frac{\partial z(\varphi ,\theta )}{r\partial \theta } =\varsigma _{\theta } (\varphi ,\theta )$$
is the slope angle. Let us introduce $x = \cos\theta$ and present the relief as a sum of normalized associated Legendre polinomials $P_{j}^{\left|m\right|} (x)$,

\begin{eqnarray} \label{13}
\frac{z(\varphi ,\theta )}{r} =-\sum _{j=1}^{\infty } \sum _{m=-j}^{j} a_{m}^{j} \exp (im\varphi )\Phi _{j}^{|m|} (x),\; \\
\Phi _{j}^{|m|} (x)=N_{j}^{|m|} P_{j}^{|m|} (x),  \nonumber
\end{eqnarray}

\begin{eqnarray*}
\int _{-1}^{1} \Phi _{j}^{|m|} (x)\Phi _{i}^{|m|} (x)dx=\delta _{ji} ,\;  N_{j}^{|m|}= \\
\left(\frac{2j+1}{2} \frac{(j-|m|)!}{(j+|m|)!} \right)^{1/2} .
\end{eqnarray*}

Then due to eq. \eqref{13} the relief slope is expressed through the derivative,

\begin{equation} \label{14}
\zeta _{\theta } (\phi ,\theta )=\sum _{j=1}^{\infty } \sum _{m=-j}^{j} a_{m}^{j} \exp (im\varphi )\sin \theta \frac{d\Phi _{j}^{|m|} (x)}{dx} .
\end{equation}

Consider the mean normalized energy (in quad\-ra\-tic sense) of the relief slopes as follows,

\begin{eqnarray} \label{15}
 E_{\theta } =\int _{0}^{2\pi } \int _{-1}^{1} < \zeta _{\theta }^{2} > d\varphi dx =
\end{eqnarray}
\begin{eqnarray}
 2\pi \sum _{j=1}^{\infty } \sum _{m=-j}^{j} <|a_{m}^{j} |^{2} >\int _{-1}^{1} (1-x^{2} )\left(\frac{d\Phi _{j}^{|m|} (x)}{dx} \right)^{2} dx.  \nonumber
\end{eqnarray}

We consider the mean values of the spectral components $<|a_{m}^{j} |^{2} >$ to be independent of $m$, which is confirmed by observations and on practice [\itc{Kaula,} \reflink{kaul66}{1966}; \itc{Rexer and Hirt,} \reflink{rexe15}{2015}; \itc{Turcotte,} \reflink{turc97}{1997}]. It thus depends only meridional index $j$ and can be found on summing over longitudinal index $m$:

\begin{equation} \label{16}
<|a_{m}^{j} |^{2} > = \frac{\alpha ^{2} }{j(j+1)(j+1/2)} ,
\end{equation}

\noindent where $\alpha$ is the mean angle of the relief slope. Taking this into account and calculating the integrals in eq. \eqref{15} we have on summing over $m$

\begin{equation} \label{17}
E_{\theta } = 2\pi \sum _{j=1}^{\infty }S_{j}  , \quad S_{j} = \frac{\alpha ^{2} }{j(j+1)} ,
\end{equation}

We thus obtain the constant spectrum of ``white\-noise'', i.e. $\delta$-correlation of $\theta$ Then for the relief energy (as mentioned above), we derive a quadratic value integrated over the spectrum from \eqref{13},

\begin{eqnarray} \label{18}
E = \int_{0}^{2\pi}\int_{-1}^{1} < z^{2} > d\varphi dx = \\
2\pi r^{2}\sum_{j-1}^{\infty}\sum_{m=-j}^{j} < |\alpha_{m}^{j}|^{2} > \int_{-1}^{1} (\Phi_{j}^{|m|}(x))^{2}dx = \nonumber  \\
2\pi r^{2}\sum_{j-1}^{\infty}\sum_{m=-j}^{j} < |\alpha_{m}^{j}|^{2} > . \nonumber
\end{eqnarray}

It is clear from here that the mean of the spectral components \eqref{17} provides the spectrum relief $\sim j^{-2}$ at $j>>1$, as was found from the empirical data of refs [\itc{Kaula,} \reflink{kaul66}{1966}; \itc{Konopliv et al.,} \reflink{kono14}{2014}; \itc{McMahon et al.,} \reflink{mcma16}{2016}; \itc{Rexer and Hirt,} \reflink{rexe15}{2015}; \itc{Turcotte,} \reflink{turc97}{1997}].

\begin{equation}\label{19}
E=2\pi \sum _{j=1}^{\infty } \sum _{m=-j}^{j} \frac{r^{2} a^{2} }{j(j+1)(j+1/2)} =  \end{equation}
\begin{eqnarray*}
2\pi \sum _{j=1}^{\infty } \frac{r^{2} a^{2} (2j+1)}{j(j+1)(j+1/2)} = \\
4\pi \sum _{j=1}^{\infty } \frac{r^{2} a^{2} }{j(j+1)} = 4\pi r^{2} \alpha ^{2} =4\pi rD_{1} ,
\end{eqnarray*}

Here $D_{1}$ is the horizontal diffusion coefficient and $\displaystyle{\sum _{j=1}^{\infty }\frac{1}{j(j+1)} =1}$ (note $\frac{1}{j(j+1)} = \displaystyle{\frac{1}{j} -\frac{1}{j+1}} $ identity). The dependence of our harmonic coefficients $1/j(j+1)$ (see eq. \eqref{19}) differs from Kaula's $1/j^{2}$ by the factor $1/j^{2}(j+1)$. At $j$ = 4 it is equal to 1/80.

In order to replace integers $j$ by continuous wave numbers $j$ at $j\to \infty $, i.e. for small horizontal scales, we shall use in eq. \eqref{13} the asymptotics of the polynoms $P_{j}^{\left|m\right|} (\cos \theta )$ at $j\to \infty :$

\begin{equation} \label{20}
P_{j}^{\left|m\right|}(\cos\theta)\approx(-1)^{m}\left(\frac{2}{\pi j\cos\theta} \right)^{1/2}\times
\end{equation}
\begin{equation}
\times \cos\left(\alpha-\frac{\pi}{4}+\frac{m\pi}{2}\right),\quad   \alpha =\left(j+\frac{1}{2} \right)\theta  \nonumber
\end{equation}

Now instead of co-latitude $\theta$ we will use the length $y$ along the meridian assuming $\theta = 0$ at the north pole,

\begin{equation} \label{21}
y = r\theta, \quad 0\le \theta < \pi.
\end{equation}

\noindent Then we rewrite the phase $\alpha$ in eq. \eqref{20} as $\alpha = ydk(j+\frac{1}{2})$, where $dk=\frac{j}{r}$ is an increment of the wave number $k = j \cdot dk$ at $j\to \infty $. Then for very small horizontal scales one has at $dk\to \infty $,
\begin{equation}\label{22}
    \alpha = ky,  \quad    k = j \cdot dk = j/r
\end{equation}

Now let us apply the wave number $k$ in the eq. \eqref{18} for the relief energy. At $j\to \infty $ we obtain

\begin{equation} \label{23}
E = 2\pi D_{1} k^{-2} = 4\pi r\alpha ^{2} k^{-2} .
\end{equation}

Here $4\pi r\alpha ^{2} k^{-2} = S(k)$ is the spectrum with dimension of the cubic length (recall that $\alpha$ the rms slope angle). It can be seen that the value $r\alpha ^{2}$ in eq. \eqref{23} is proportional to the diffusion coefficient $D_{1}$ in the coordinate space.

The relief spectrum in eq. \eqref{19} differs from the Kaula's ``rule of the thumb'' $\propto j^{-2} $ by the multiplier $[j(j + 1)]^{-1}$. The relative difference between the two at $j\to \infty $ decreases as $[j^{2}(j + 1)]^{-1}$ but for small $j$ it is quite noticeable, e.g. at $j = 1$ the difference is 1/2, and at $j = 4$ it is 20\%. This shows that the hypothesis on the ``white noise'' of the relief slope angles is not accurate at these scales. Thus the corresponding relief structure at $j\le 5$ could be different for each celestial body due to its peculiar tectonics.



\section{Results and Discussion}

Our scheme of the relief evolution basing on eq. \eqref{9} of the probability density changes is the horizontal diffusion of the tectonically formed vertical structure under the gravity acting along slopes. Slopes as a measure of forces resist to winds; water and rocks are running down along them, etc. \figref{2} (Figure 7.19 from [\itc{Turcotte,} \reflink{turc97}{1997}]) presents the amplitudes of the relief spherical harmonics: for Earth $j\le 180$ and for Venus $j\le 50$. At $j\ge 4$ all coefficients approach to the inverse square dependence.

We compare our theoretical expression for harmonic coefficients from eq. \eqref{19}

\begin{equation} \label{24}
S_{j} =\frac{4\pi rD_{1} }{j(j+1)}
\end{equation}
to their empirical values from \figref{2}.

\begin{figure}                      % Fig 2
\figurewidth{90mm}
\setimage{}{}{}{140mm}{2019es000651-f02}
\fcaplabel{Relief spectral harmonics for Earth and Venus after [\itc{Turcotte,} \reflink{turc97}{1997}].}
\end{figure}


The comparison has been done for selected harmonics with accuracy about 10\% within their value order. The harmonic numbers for Earth were 5, 10, 20, 30, 40, 60, 90. This produced $D_{1} = 1.3 \pm 0.3$ m. For Venus the numbers 5, 10, 15, 20, 30, 40 and 50 produced $D_{1} = 0.19 \pm 0.03$ m. Our diffusion coefficient is determined with accuracy about 20\%. The striking difference between these two planets has been already noted by [\itc{Turcotte,} \reflink{turc97}{1997}]. However, one should recall that Venusian data are related only to equatorial plates [\itc{Turcotte,} \reflink{turc97}{1997}]. The global relief characteristics can be estimated with much higher precision and efforts (see e.g. [\itc{Rexer and Hirt,} \reflink{rexe15}{2015}]), which we leave for much younger colleagues, but we believe that it would be worth doing.

\figref{3} presents the relief spectra as the functions of the surface linear wave length $k$.  The authors of [\itc{Rexer and Hirt,} \reflink{rexe15}{2015}] have computed spherical harmonics for Earth and Moon up to $j$ = 46200 and for Mars up to $j$ = 23100. Our eqs. \eqref{20} -- \eqref{23} present high frequency harmonics in terms of $k^{-2} $. At \figref{3} we added straight lines corresponding to the spectra $k^{-2} $ and $k^{-4} $, the last for very high wave numbers. The dependence $k^{-2} $ holds well for Mars, line 1, and for three terrestrial spectra, lines 3-5, with some steeping for waves shorter than a few kilometers. The Moon spectrum falls down and remains almost flat at 200 -- 80 km, being very specific in this respect. This evidences that the Moon is bombarded by the asteroids the size of about 10 km, since the size of a crater is 10-20 times exceeds the diameter of the asteroid.

The explanation for the relief spectral decrease from $k^{-2} $ to $k^{-4} $ can be explained by rejecting the $\delta$-correlation for the relief slopes hypothesis. The simplest correlation function is



\begin{sidewaysfigure}                      % Fig 3
\figurewidth{140mm}
\setimage{}{}{155mm}{}{2019es000651-f03}
\fcaplabel{Relief spectra as functions of the surface wavelength for Moon, Mars and Earth [\itc{Rexer and Hirt,} \reflink{rexe15}{2015}](3 cases labeled at lower left). We added lines $k^{-2}$ and $k^{-4}$ for large and small distances.}
\end{sidewaysfigure}

\begin{equation} \label{25}
B_{\theta } (\theta (y_{1} )\theta (y_{2} ))\sim \exp (-\beta y), \quad y=y_{2} -y_{1},
\end{equation}

\noindent where $\beta \approx 1/h$ and $h$ is mean relief height. Its spectrum is

\begin{equation} \label{26}
S_{\theta } (k)=\frac{\beta }{\pi (\beta ^{2} +k^{2} )} ,
\end{equation}

\noindent and if we know the spectrum of the derivatives, i.~e. angles, the spectrum of the value itself, relief, would be

\begin{equation} \label{27}
S_{\theta } (k)=k^{2} S_{\theta }^{2} (k)=\frac{k^{-2} \beta }{\pi (\beta ^{2} + k^{2} )} .
\end{equation}

This spectrum has both right asymptotes: at $h^{2} \ll \beta ^{2} $ it would be proportional to $k^{-2} $ and in inverse case for small distances it $\sim k^{-4} $. The gravity acceleration on the Moon is six times weaken that on Earth and other conditions equal the relief inhomogeneities there could be formed easier and higher.



\section{Conclusion}

Our analysis is based on the Fokker -- Planck -- Kolmogorov, FPK, equation cast into eq. \eqref{9} for altimeter records. It has the solution $\left\langle h^{2} \right\rangle = D_{1} y$ for the mean square of the relief height. Here $D_{1}$ is the horizontal diffusivity of the relief. Such a solution has the spectrum $S(k)=Dk^{-2}$, with $k = 2\pi /\lambda$ along the horizontal direction. The modified FPK equation is based on the hypothesis of the time $\delta$-correlation of acting forces which in this case corresponds to horizontal $\delta$-correlation of the relief slope angles. The spherical harmonics of the relief have been found by Kaula to follow the inverse-square rule for the harmonic number $j$. Our ``white noise'' hypothesis on the horizontal spectrum of the slopes slightly modifies the Kaula's rule for high harmonics by replacing the inverse square dependence by the factor $[j(j + 1)]^{-1}$. The difference decreases as $[j^{2}(j + 1)]^{-1}$. The numerator in eq. \eqref{24} gives the constant for the Kaula's rule.

Rather crude estimate of diffusion coefficients $D_{1}$ performed for Earth yields $D =~1.3~$m, several times higher than for Venus. A similar investigation may also be performed for the planetary gravity field fluctuations with similar results, though in this case one should account for the internal structure and tectonics. Consideration of statistical scales \eqref{4} as structure functions with zero initial conditions allows us to propose their spectral representation. For the very high harmonics with continuous spectrum the spectral behaviour changes from $k^{-2} $ to $k^{-4} $, see eqs. \eqref{25} -- \eqref{27}.

It can be concluded that the Kaula's rule for gravity field as well as for the relief is an asymptotic consequence of the random walk laws, for the first time formulated by Einstein (see [\itc{Bridgman,} \reflink{brid21}{1921}]) and in the most general and practical way by Kolmogorov and his successors in 1934 [\itc{Monin and Yaglom,} \reflink{moni75}{1975}; \itc{Obukhov,} \reflink{obuk59}{1959}; \itc{Yaglom,} \reflink{yagl86}{1986}].

Thus the line of study may be considered as the model, the base for use the dimensional analysis for investigation of various processes. When obtaining dimensional formulas one always should remember that the formulas must contain a numerical coefficient which could be retrieved from the data of observations. The formulas can be valid only in a certain range of dimensionless similar numbers. In certain situations there can be more than one such number. How should one act in these cases? One should consult to professor \itc{G. I. Barenblatt,} \reflink{bare96}{1996}; \reflink{bare03}{2003}], who was Kolmogorov's student.

Geophysical interpretation of these results is waiting for a thorough analysis. Because the Kaula's constant for gravity field is used to obtain the information on the internal structure of celestial bodies (mascones, isostasy), the diffusion coefficient $D_{1}$ could contain the information on tectonics, surface material properties, etc. This coefficient is an analogue for the Kaula's gravity field fluctuation constant. Much work has yet to be done by much younger people since the total age of the two present authors is not far from 160 yrs.

\acknowledgment{We thank A.~A. Lushnikov for carefully reading the manuscript and indicating some errors in the formulas.  This work has been particularly supported by the RAS Presidium Program No. 7 ``The development of non-linear methods in theoretical and mathematical physics''. The authors are very thankful to Dr. O. G. Chkhetiany for discussion of the results and his continuous help during the long work.}

\begin{references}
\reference{bare96} Barenblatt, G. I. (1996) , {\itshape Scaling, Similarity and Intermediate Asymptotics}, 386 pp., Cambridge Univ. Press, Cambridge, UK,  \href{https://doi.org/10.1017/CBO9781107050242}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{bare03} Barenblatt, G. I. (2003) , {\itshape Scaling}, 171 pp., Cambridge Univ. Press, Cambridge, UK,  \href{https://doi.org/10.1017/CBO9780511814921}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{brid32} Bridgman, P. W. (1932) , {\itshape Dimensional Analysis}, 113 pp., Yale Univ. Press, London, UK.

\reference{gled10} Gledzer, E. B., G. S. Golitsyn (2010) , Scaling and finite ensembles of particles with the energy influx, {\itshape Doklady RAS}, {\itshape 433}, no. 4, p. 466--470.

\reference{goli01} Golitsyn, G. S. (2001) , The place of the Gutenberg-Richter law among other statistical laws of nature, {\itshape Comput. Seismology}, {\itshape 77}, p. 160--192.

\reference{goli03} Golitsyn, G. S. (2003) , Statistical description of the planetary surface relief and its evolution, {\itshape Proc. RAS. Physics of the Earth}, no. 7, p. 3--8.

\reference{goli18} Golitsyn, G. S. (2018) , Laws of random walks by A. N. Kolmogorov 1934, {\itshape Soviet Meteorol. Hydrol}, no. 3, p. 5--15,  \href{https://doi.org/10.3103/S1068373918030019}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{goli18} Golitsyn, G. S. (2018) , Random Walk Laws by A.~N. Kolmogorov as the bases for understanding most phenomena of the nature, {\itshape Izv. RAS Atmos. and Ocean. Phys.}, {\itshape 54}, p. 223--228.

\reference{kaul66} Kaula, W. M. (1966) , {\itshape Theory of Satellite Geodesy}, 143 pp., Bleinsdell, Waltham. Ma.

\reference{kolm34} Kolmogorov, A. N. (1934) , Zufallige Bewegungen, {\itshape Ann. Math. }, {\itshape 35}, p. 166--117,  \href{https://doi.org/10.2307/1968123}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{kono14} Konopliv, A. S., et al. (2014) , The Vesta gravity field, {\itshape Icarus}, {\itshape 240}, p. 103--117,  \href{https://doi.org/10.1016/j.icarus.2013.09.005}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{lind08} Lindborg, E. (2008) , Stratified turbulence: a possible interpretation of some geophysical turbulence measurements, {\itshape J. Atmos. Sci.}, {\itshape 65}, p. 2416--2424,  \href{https://doi.org/10.1175/2007JAS2455.1}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{mcma16} McMahon, J. W., et al. (2016) , Understanding the Kaula's rule for small bodies, Proc. 47th Lunar and Planetary Sci. Conference, AGU, Washington, DC (\href{https://www.hou.usra.edu/meetings/lpsc2016/pdf/2129.pdf}{https://www.hou.usra.edu/\linebreak meetings/lpsc2016/pdf[linebreak] 2129.pdf}).

\reference{moni75} Monin, A. S., A. M. Yaglom (1975) , {\itshape The Statistical Hydromechanics}, MIT Press, Cambridge (Russ. ed. 1967).

\reference{obuk41} Obukhov, A. M. (1941) , On the energy distribution in the spectrum of turbulent flow, {\itshape Proc. (Doklady) USSR Ac. Sci.}, {\itshape 32}, no. 1, p. 22--24.

\reference{obuk59} Obukhov, A. M. (1959) , Description of turbulence in terms of Lagrangian variables, {\itshape Adv. Geophys.}, {\itshape 6}, p. 113--115,  \href{https://doi.org/10.1016/S0065-2687(08)60098-9}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{rexe15} Rexer, M., C. Hirt (2015) , Ultra-high spherical harmonics analysis and application to planetary topography of Earth, Mars and Moon, {\itshape Surv. Geophys.}, {\itshape 36}, no. 6, p. 803--830,  \href{https://doi.org/10.1007/s10712-015-9345-z}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{rich26} Richardson, L. F. (1926) , Atmospheric diffusion shown on a distance-neighbour graph, {\itshape Proc. Roy. Soc.}, {\itshape A110}, no. 756, p. 709--737,  \href{https://doi.org/10.1098/rspa.1926.0043}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{rich29} Richardson, L. F. (1929) , A search for the law of atmospheric diffusion, {\itshape Beitr. Phys. f\"{u}r Atmosph\"{a}re}, {\itshape 15}, no. 1, p. 24--29.

\reference{turc97} Turcotte, D. L. (1997) , {\itshape Chaos and Fractals in Geology and Geophysics}, 416 pp., Cambridge University Press, New York,  \href{https://doi.org/10.1017/CBO9781139174695}{\textcolor{brown}{\space\bfseries Cross\-ref}}

\reference{yagl55} Yaglom, A. M. (1955) , Correlation theory for processes with stochastic stationary increments of the n-th order, {\itshape Math. Collection}, {\itshape 37}, no. 1, p. 259--288 (in Russian).

\reference{yagl86} Yaglom, A. M. (1986) , {\itshape Correlation Theory of Stationary and Related Random Functions}, 526 pp., Springer Verlag, Berlin.

\end{references}
\end{document}
