text:
# Correlations and screening of topological charges in gaussian random fields

## 1 Introduction

Although a spatially extended field may be smooth, and contain no infinities, it may nonetheless have point singularities associated with its topology. For example, a 2-dimensional landscape, specified by a real function, has *critical points* (stationary points), where the gradient of the field vanishes and the gradient direction cannot be defined; in the neighbourhood of such a point, the gradient direction changes arbitrarily quickly through all of its values. Such points are characterised by a *topological charge*, a signed integer which is determined by the local geometry of the singularity; in the case of critical points, the number is the (signed) number of rotations of the gradient vector around the critical point; it is $`+1`$ for maxima and minima, and $`1`$ for saddles. The topological charge is also called the *Poincaré index* (or *Hopf index*) (Milnor 1965), and is defined as a signed integer at the point zeros of $`n`$-dimensional real vector fields in $`n`$ dimensions. Although the number of zeros may change as a field evolves, the total topological charge is constant; it is a topological invariant.

Here, I discuss the statistical properties of various zeros/singularities, in fields which are specified by gaussian random functions. In this case, only structurally stable zeros of charge $`\pm 1`$ occur. Specifically, I shall discuss how densities and two point charge correlation functions of the distributions of signed singularities may be calculated under a rather general scheme, and then use this scheme to calculate the topological charge correlation functions of three types of singularity in isotropic random fields: zeros of $`n`$-dimensional vectors in $`n`$-dimensional space (rederiving a result originally due to Halperin (1981)); critical points of random scalar fields in two dimensions; and umbilic points of random scalar fields in two dimensions. The charge correlation functions for each of these types of topological singularity are found to be different in each case, and dependent on the underlying correlation function of the gaussian field. Each function is found to satisfy a screening relation associated with ionic liquids.

Topological zeros are very important in many areas of physics and mathematics: in addition to critical points (gradient zeros) which have obvious importance, the zeros of 2-dimensional complex scalar fields (phase singularities, wave dislocations, vortices, realised as 2-dimensional vectors) are also of great interest, especially where the field is a quantum wavefunction or an optical field (e.g. Nye and Berry 1974, Berry and Dennis 2000, Dennis 2001b), or an order parameter (Mermin 1979). Only point zeros are considered here - the dimensionality $`n`$ of the vector field, whose zeros provide the singularity, is assumed equal to the dimensionality of the space.

The first systematic study of the statistical geometry of random real scalar fields in two dimensions was by Longuet-Higgins (1957a,b, 1958), who generalised one-dimensional methods of Rice (1954) to calculate, amongst other things, the densities of critical points and probability density function of the gaussian curvature of the function. Halperin (1981) derived the $`n`$-dimensional vector correlation function, whose proof was supplied by Liu and Mazenko (1992), and has recently been recast in the language of riemannian geometry by Foltin (2003a). Various correlation functions in two dimensions, including those for phase singularities and critical points, were investigated numerically by Freund and Wilkinson (1998). Planar phase singularity correlations (including density correlations) were investigated by Berry and Dennis (2000).

The topological singularities of interest are the zeros of an $`n`$-dimensional vector field $`\mathrm{v}\mathbf{=}\mathbf{(}\mathrm{v}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_\mathrm{n}\mathbf{)}\mathbf{.}`$ The field is a smooth function defined on an $`n`$-dimensional euclidean space, with points labelled by the vector $`\mathrm{r}\mathbf{=}\mathbf{(}\mathrm{r}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{r}_\mathrm{n}\mathbf{)}\mathbf{.}`$ In random fields, the only statistically significant zeros are those of first order, whose jacobian determinant

$$𝒥=det_iv_j,$$
(1.1)
is nonzero (where $`_i_{,i}/r_i`$). The topological charge of such zeros is given by $`\mathrm{sign}𝒥.`$

When the field is random, the density $`d(\mathrm{r}\mathbf{)}`$ of zeros at a position $`\mathrm{r}`$ is

$$d(\mathrm{r}\mathbf{)}\mathbf{=}\mathbf{}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}\mathbf{)}\mathbf{)}\mathbf{|}𝓙\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}\mathbf{)}\mathbf{)}\mathbf{|}\mathbf{}\mathbf{,}$$
(1.2)
where $``$ denotes averaging over the statistical ensemble, and the zeros are identified by the $`n`$-dimensional Dirac $`\delta `$-function. The modulus of the jacobian is included to ensure that each zero has the correct statistical weight, so $`d(\mathrm{r}\mathbf{)}`$ has the units of density. This expression confirms that degenerate zeros are so rare they have no statistical significance. The average topological charge $`q(\mathrm{r}\mathbf{)}`$ at $`\mathrm{r}`$ is expressed as the density (1.2), but each singularity is weighted by its charge, that is, the sign of the jacobian. The charge correlation functions $`g`$ with which this paper is concerned are the generalisation of the average charge density to that at two points $`\mathrm{r}_\mathrm{A}\mathbf{,}\mathrm{r}_\mathrm{B}\mathbf{,}`$ normalised by the density:

$$g(\mathrm{r}_\mathrm{A}\mathbf{,}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{=}\frac{1}{\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}}\mathbf{}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}_\mathrm{A}\mathbf{)}\mathbf{)}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{)}𝓙\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}_\mathrm{A}\mathbf{)}\mathbf{)}𝓙\mathbf{(}\mathrm{v}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{)}\mathbf{}\mathbf{.}$$
(1.3)
In the following section, the scheme for calculating topological charge correlation functions (1.3) in general gaussian random fields is explained. Section 3 then provides background to gaussian fields which are statistically stationary and isotropic. The scheme is then applied to zeros of $`n`$-dimensional vector fields (section 4), critical points of 2-dimensional scalar fields (section 5), and umbilic points of the same fields (section 6). The phenomenon of topological charge screening in these three cases is discussed in section 7.

As this paper was being completed, I became aware of the work of Foltin (2003b), who performs similar calculations for critical points by a different, possibly simpler method.

## 2 Gaussian random functions and scheme for calculating charge correlation functions

An ensemble of scalar functions $`f`$ in $`n`$-dimensional space is said to be a gaussian random function if the probability distribution of the function at each point of space is given by a gaussian distribution, so the function at each point defines a gaussian random variable (Adler 1981). The only restrictions on gaussian random functions made in this section are that they be centred (the average $`f=0`$) and that first derivatives exist.

Our starting point is the well-known expression for the probability density function for a set of $`N`$ independent gaussian random functions $`\mathrm{u}\mathbf{=}\mathbf{(}\mathrm{u}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{u}_\mathrm{N}\mathbf{)}\mathbf{,}`$

$$P(\mathrm{u}\mathbf{)}\mathbf{=}\frac{\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{u}\mathbf{}\mathrm{\Sigma }^\mathbf{}1\mathbf{}\mathrm{u}\mathbf{/}2\mathbf{)}}{\mathbf{(}2𝝅\mathbf{)}^{\mathrm{N}\mathbf{/}2}\sqrt{𝐝𝐞𝐭\mathrm{\Sigma }}}\mathbf{,}$$
(2.1)
where $`𝚺`$ is the correlation matrix, with components defined

$$\mathrm{\Sigma }_{ij}=u_iu_j.$$
(2.2)
It is now possible to present the general scheme for calculating 2-point topological charge correlation functions of zeros of gaussian random functions, expressed in (1.3). The topological charges considered are all zeros of a real vector gaussian random function $`\mathrm{v}\mathbf{=}\mathbf{(}\mathrm{v}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_\mathrm{n}\mathbf{)}\mathbf{,}`$ of dimension $`n,`$ and the correlation function is the average of the product of the local density at points labelled $`A,B.`$ Quantities evaluated at these places are denoted with the appropriate subscript $`A`$ or $`B.`$ The variables that appear in the average (1.3) are therefore the components of $`\mathrm{v}_\mathrm{A}\mathbf{,}\mathrm{v}_\mathrm{B}\mathbf{,}`$ and their derivatives $`\mathrm{v}_\mathrm{A}\mathbf{,}\mathbf{}\mathrm{v}_\mathrm{B}`$ that appear in the jacobians $`𝒥(\mathrm{v}_\mathrm{A}\mathbf{)}\mathbf{}𝓙_\mathrm{A}\mathbf{,}𝓙\mathbf{(}\mathrm{v}_\mathrm{B}\mathbf{)}\mathbf{}𝓙_\mathrm{B}\mathbf{.}`$ There are a total of $`m`$ different derivatives $`_iv_j`$ appearing in each jacobian, and in general $`mn^2.`$ For example, in the case of critical points, $`n=2`$, and $`m=3`$ (the terms appearing in the jacobian are $`_{11}f,_{22}f,_{12}f=_{21}f`$). The calculation therefore involves an average of the $`N=2(m+n)`$-dimensional gaussian random vector

$$\mathrm{u}\mathbf{=}\mathbf{(}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{A}\mathbf{)}_1\mathbf{,}\mathbf{}\mathbf{,}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{A}\mathbf{)}_\mathrm{m}\mathbf{,}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{B}\mathbf{)}_1\mathbf{,}\mathbf{}\mathbf{,}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{B}\mathbf{)}_\mathrm{m}\mathbf{,}\mathrm{v}_{\mathrm{A}1}\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_{\mathrm{A}\mathrm{n}}\mathbf{,}\mathrm{v}_{\mathrm{B}1}\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_{\mathrm{B}\mathrm{n}}\mathbf{)}\mathbf{.}$$
(2.3)
The correlation matrix $`𝚺`$ for this $`\mathrm{u}`$ is defined as in (2.2), and averages are evaluated according to the probability density function (2.1).

The correlation function in (1.3) is therefore

$`g_{AB}={\displaystyle \frac{\delta ^n(\mathrm{v}_\mathrm{A}\mathbf{)}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}_\mathrm{B}\mathbf{)}𝓙_\mathrm{A}𝓙_\mathrm{B}\mathbf{}}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}}}`$
$`={\displaystyle \frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^{\mathrm{n}\mathbf{+}\mathrm{m}}\sqrt{𝐝𝐞𝐭\mathrm{\Sigma }}}}`$
$`\times {\displaystyle }^{2(n+m)}\mathrm{u}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}_\mathrm{A}\mathbf{)}𝜹^\mathrm{n}\mathbf{(}\mathrm{v}_\mathrm{B}\mathbf{)}𝓙_\mathrm{A}𝓙_\mathrm{B}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{u}\mathbf{}\mathrm{\Sigma }^\mathbf{}1\mathbf{}\mathrm{u}\mathbf{/}2\mathbf{)}\mathbf{.}`$ (2.4)
where $`d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathbf{,}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}`$ are the appropriate densities of zeros (1.2). The jacobians $`𝒥_A,𝒥_B`$ can be quite complicated; each is a multilinear function, involving a sum of $`n`$-fold products of $`u_1,\mathrm{},u_m`$ (for $`A`$) or $`u_{m+1},\mathrm{},u_{2m}`$ (for $`B`$).

The $`2n\times 2n`$-dimensional lower right submatrix of $`𝚺`$ (i.e. the averages dependent only on $`\mathrm{v}_\mathrm{A}\mathbf{,}\mathrm{v}_\mathrm{B}`$) is denoted by $`𝐊;`$ the complementary $`2m\times 2m`$-dimensional upper left submatrix of $`𝚺^1`$ is denoted by $`𝚵^1;`$ that is,

$$𝚺=\left(\begin{array}{cc}& \\ & \\ & 𝐊\end{array}\right),𝚺^1=\left(\begin{array}{cc}𝚵^1& \\ & \\ & \end{array}\right),$$
(2.5)
where the other blocks marked $``$ are not denoted by special symbols. In A, Jacobi’s determinant theorem is used to show that

$$det𝚺=det𝐊det𝚵.$$
(2.6)
The integrals in $`v_{A1},\mathrm{},v_{Bn},`$ only involving $`\delta `$-functions, are performed, leaving the $`2m`$ derivative terms $`\mathrm{v}`$ to integrate; let $`\mathrm{u}^{\mathbf{}}\mathbf{=}\mathbf{(}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{A}\mathbf{)}_1\mathbf{,}\mathbf{}\mathbf{,}\mathbf{(}\mathbf{}\mathrm{v}_\mathrm{B}\mathbf{)}_\mathrm{m}\mathbf{)}\mathbf{.}`$ Then

$`g_{AB}={\displaystyle \frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^{\mathrm{n}\mathbf{+}\mathrm{m}}\sqrt{𝐝𝐞𝐭\mathrm{\Sigma }}}}{\displaystyle ^{2m}\mathrm{u}^{\mathbf{}}𝓙_\mathrm{A}𝓙_\mathrm{B}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{u}^{\mathbf{}}\mathbf{}\mathrm{\Xi }^\mathbf{}1\mathbf{}\mathrm{u}^{\mathbf{}}\mathbf{/}2\mathbf{)}}`$
$`={\displaystyle \frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^{\mathrm{n}\mathbf{+}\mathrm{m}}\sqrt{𝐝𝐞𝐭\mathrm{\Sigma }}}}`$
$`\times {\displaystyle }^{2m}\mathrm{u}^{\mathbf{}}𝓙_\mathrm{A}𝓙_\mathrm{B}{\displaystyle \frac{\sqrt{𝐝𝐞𝐭\mathrm{\Xi }}}{\mathbf{(}2𝝅\mathbf{)}^\mathrm{m}}}{\displaystyle \mathbf{}}\mathbf{}^{2\mathrm{m}}\mathrm{t}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{t}\mathbf{}\mathrm{u}^{\mathbf{}}\mathbf{}\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}\mathbf{,}`$ (2.7)
upon Fourier transforming the gaussian in $`\mathrm{u}^{\mathbf{}}`$ with the $`2m`$-dimensional Fourier vector variable $`\mathrm{t}\mathbf{.}`$ The jacobian terms $`𝒥_A`$ (depending on $`u_1^{},\mathrm{},u_m^{}`$) and $`𝒥_B`$ (depending on $`u_{m+1}^{},\mathrm{},u_{2m}^{}`$) may be replaced with partial derivative terms in $`t_j,`$

$$u_j^{}\mathrm{}_{t_j}=\mathrm{}_j,$$
(2.8)
where $``$ is used to denote partial derivatives in $`2m`$-dimensional $`\mathrm{t}`$-space. The jacobians, with this replacement, have now become differential operators, denoted $`𝒥_A,𝒥_B.`$ Therefore

$`g_{AB}={\displaystyle \frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^{\mathrm{n}\mathbf{+}2\mathrm{m}}}}\sqrt{{\displaystyle \frac{det𝚵}{det𝚺}}}`$
$`\times {\displaystyle }^{2m}\mathrm{t}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}𝓙_\mathbf{}\mathrm{A}𝓙_\mathbf{}\mathrm{B}{\displaystyle \mathbf{}}\mathbf{}^{2\mathrm{m}}\mathrm{u}^{\mathbf{}}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{t}\mathbf{}\mathrm{u}^{\mathbf{}}\mathbf{)}`$
$`={\displaystyle \frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^{\mathrm{n}\mathbf{+}2\mathrm{m}}}}{\displaystyle \frac{(2\pi )^{2m}}{\sqrt{det𝐊}}}`$
$`\times {\displaystyle }^{2m}\mathrm{t}\mathrm{𝐞𝐱𝐩}\mathbf{(}\mathbf{}\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}𝓙_\mathbf{}\mathrm{A}𝓙_\mathbf{}\mathrm{B}𝜹^{2\mathrm{m}}\mathbf{(}\mathrm{t}\mathbf{)}\mathbf{,}`$ (2.9)
where, in the second line, the integral over $`\mathrm{u}^{\mathbf{}}`$ is realised as the Fourier transform of the $`\delta `$-function of $`\mathrm{t}\mathbf{,}`$ and (2.6) has been used in the prefactor. The expression is then integrated by parts (so the $`𝒥_{}`$ operators act on the exponential term rather than the $`\delta `$-function), and then integrated in $`\mathrm{t}\mathbf{.}`$ The final expression is

$$g_{AB}=\frac{1}{d(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{d}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{(}2𝝅\mathbf{)}^\mathrm{n}}\frac{1}{\sqrt{det𝐊}}D,$$
(2.10)
where

$$D=\left[𝒥_A𝒥_B\mathrm{exp}(\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}\right]_{\mathrm{t}\mathbf{=}\mathrm{𝟎}}.$$
(2.11)
The charge correlation function therefore only depends on $`d,det𝐊`$ and components of the inverse reduced inverse correlation matrix $`𝚵.`$

This final step, evaluating $`D`$ in (2.11) is the most complicated part of the calculation, and depends on the precise form of the jacobian determinant $`𝒥.`$ Each summand in the operator $`𝒥_A𝒥_B`$ is a $`2n`$-fold derivative over some index set $`\{\tau \}=\{\tau _1,\mathrm{},\tau _{2n}\},`$ and is comprised of $`n`$ terms from $`_1,\mathrm{},_m`$ (for $`A`$), and $`n`$ terms from $`_{m+1},\mathrm{},_{2m}`$ (for $`B`$), either set possibly including repetitions. The result of this particular operation is

$$\left[\left(\underset{j=1}{\overset{2n}{}}\mathrm{}_{\tau _j}\right)\mathrm{exp}(\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}\right]_{\mathrm{t}\mathbf{=}\mathrm{𝟎}}=\underset{p}{\overset{\mathrm{all}\mathrm{pairings}}{}}\underset{(\tau _\lambda ,\tau _\lambda ^{})p}{}\mathrm{\Xi }_{\tau _\lambda \tau _\lambda ^{}},$$
(2.12)
where the sum on the right hand side is over all pairings $`p=(\tau _\alpha ,\tau _\alpha ^{})\mathrm{}(\tau _\mu ,\tau _\mu ^{})`$ of indices in $`\{\tau \};`$ there are $`(2n1)!!`$ such pairings, each involving $`n`$ sets of pairs. The $`n`$-fold product on the right hand side is over all components of $`𝚵`$ with indices given by the appropriate pair. The $`\mathrm{\Xi }_{ij}`$ themselves are found by further application of Jacobi’s determinant theorem in A, and are expressed in terms of minors of $`det𝚺`$ in (1.5).

The correlation functions calculated explicitly in this paper are not too complicated, either because $`n`$ is small (only two-dimensional fields are considered in sections 5, 6), or $`𝚵`$ is sparse, as in the case of random vectors (section 4).

The density $`d`$ of zeros appears in (2.10); in general, this may be difficult to calculate, due to the modulus sign in (1.2). For random $`n`$-dimensional vector fields, the main part of the density calculation is in B. For critical and umbilic points, these densities were calculated by Longuet-Higgins (1957a,b) and Berry and Hannay (1977).

## 3 Isotropic random fields

The topological charge correlation function in (2.10) is extremely general, applying to any centred differentiable gaussian random vector field. In this section, and for the remainder of the article, attention will be restricted to stationary isotropic random fields. For these fields, all averages are (statistically) translation and rotation invariant. They are conveniently given by a Fourier representation

$$f(\mathrm{r}\mathbf{)}\mathbf{=}\underset{\mathrm{k}}{\mathbf{}}\mathrm{a}\mathbf{(}\mathrm{k}\mathbf{)}\mathrm{𝐜𝐨𝐬}\mathbf{(}\mathrm{k}\mathbf{}\mathrm{r}\mathbf{+}\mathit{\varphi }_\mathrm{k}\mathbf{)}\mathbf{,}$$
(3.1)
where $`\mathrm{k}`$ are now the Fourier variable vectors (wavevectors). The components $`v_i`$ of random vectors are specified by independent identically distributed realisations of (3.1). The amplitude $`a(k)`$ only depends on the magnitude $`|\mathrm{k}\mathbf{|}\mathbf{=}\mathrm{k}\mathbf{,}\mathrm{s}`$ and the phase $`\varphi _\mathrm{k}`$ is uniformly random - each ensemble member is therefore labelled by the choice of $`\varphi _\mathrm{k}`$ for each $`\mathrm{k}\mathbf{.}`$ It also may represent the spatial part of a real linear homogeneous nondispersive wavefield, for which the representation (3.1) is particularly evocative. The infinite $`\mathrm{k}`$ set is assumed sufficiently dense that they may be represented as an integral, and

$$\underset{\mathrm{k}}{}a^2(k)^n\mathrm{k}\mathrm{\Pi }\mathbf{(}\mathrm{k}\mathbf{)}\mathbf{}\mathbf{,}$$
(3.2)
where $`\mathrm{\Pi }(k)`$ is the power spectrum of the field; by the Wiener-Khinchine theorem (Feller 1950), $`\mathrm{\Pi }(k)`$ is the $`n`$-dimensional Fourier transform of the *field correlation function* $`C(r),`$ where $`r=|\mathrm{r}_\mathrm{A}\mathbf{}\mathrm{r}_\mathrm{B}\mathbf{|}\mathbf{,}`$ and

$$C(r)=f_Af_B=f(\mathrm{r}_\mathrm{A}\mathbf{)}\mathrm{f}\mathbf{(}\mathrm{r}_\mathrm{B}\mathbf{)}\mathbf{}\mathbf{,}$$
(3.3)
normalised such that $`C(0)=f^2=1.`$ The only condition on $`C`$ is that it is symmetric and has positive Fourier transform.

Averages of derivatives of $`f`$ may be represented as moments of $`\mathrm{\Pi }`$ (as by Longuet-Higgins 1957a,b, Berry and Hannay 1977, Berry and Dennis 2000), or equivalently derivatives of $`C`$ as is done here. Coordinates are chosen where

$$\mathrm{r}\mathbf{=}\mathrm{r}_\mathrm{B}\mathbf{}\mathrm{r}_\mathrm{A}\mathbf{,}\mathrm{r}_1\mathbf{=}\mathrm{r}\mathbf{,}\mathrm{r}_\mathrm{j}\mathbf{=}0\mathbf{,}\mathrm{j}\mathbf{=}2\mathbf{,}\mathbf{}\mathbf{,}\mathrm{n}\mathbf{.}$$
(3.4)
Since the fields are isotropic, the results are not affected by this choice. The correlations computed in sections 5, 6 are in two dimensions; in this case, direction 1 is denoted by $`x,`$ direction 2 by $`y.`$

Representing derivatives by subscripts, the correlations of first derivatives of $`f`$ are found to be (Berry and Dennis 2000)

$`E`$ $`f_{A,1}f_B=f_Af_{B,1}`$ $`=\left[_1C\right]_{r_1=r,r_2,\mathrm{}=0}=C^{}`$
$`F`$ $`f_{A,1}f_{B,1}=f_{A,11}f_B=f_Af_{B,11}`$ $`=\left[_1^2C\right]_{r_1=r,r_2,\mathrm{}=0}=C^{\prime \prime },`$
$`H`$ $`f_{A,j}f_{B,j}=f_{A,jj}f_B=f_Af_{B,jj}`$ $`=\left[_2^2C\right]_{r_1=r,r_2,\mathrm{}=0}=C^{}/r,j1.`$ (3.5)
Averages involving $`f`$ and its first derivatives other than those in (3.3), (3.5) are zero. The averages are equal to appropriate derivatives of $`C(\mathrm{r}\mathbf{)}\mathbf{,}`$ and then setting $`r_1=r,r_2,\mathrm{}=0`$ (as in (3.4)), and derivatives in $`r_j`$ of odd order vanish. The functions in (3.5), when $`\mathrm{r}\mathbf{=}\mathrm{𝟎}`$ (denoted by subscript 0), are

$$E_0=0,F_0=H_0=C^{\prime \prime }(0)=C_0^{\prime \prime }.$$
(3.6)
It is easily verified that $`C_0^{\prime \prime }>0,`$ since $`C`$ has a positive Fourier transform. The correlation matrices (5.4), (6.5) required to calculate the topological charge correlation functions of critical points and umbilic points, involve higher derivatives of $`C,`$ given in (5.5), (6.6).

As the separation between $`A`$ and $`B`$ increases, the correlation between $`f_A`$ and $`f_B`$ decreases and $`C0.`$ This decay is slowest in the case where all the $`\mathrm{k}`$ in (3.1) have the same length $`k_0,`$ (set to 1 for convenience) and the power spectrum $`\mathrm{\Pi }(k)=\delta (k1).`$ For any $`n>2,`$ the corresponding correlation function is a Bessel function times a factor of $`r,`$ with oscillatory decay that falls off like $`r^{(n1)/2}.`$ Of particular interest is the $`n=2`$ case, for which $`C(r)=J_0(r).`$ The spectrum in this case was called the ring spectrum by Longuet-Higgins (1957a,b) and Berry and Dennis (2000), and is conjectured to model the high eigenfunctions in quantum chaotic billiards (Berry 1978, 2002).

## 4 Correlations of zeros of isotropic vector fields

In the present section, we shall consider $`n`$-dimensional gaussian random vector fields $`\mathrm{v}\mathbf{=}\mathbf{(}\mathrm{v}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_\mathrm{n}\mathbf{)}`$ in $`n`$ dimensions whose cartesian components are independent and identically distributed gaussian random fields (the derivatives of the components are also assumed completely independent). The jacobian $`𝒥,`$ whose sign determines the topological charges of the zeros, is the determinant of the matrix of first derivatives (1.1).

We begin by calculating the density (1.2) of zeros of random vectors. This was calculated for $`n=1,2,3`$ by Halperin (1981) and Liu and Mazenko (1992), (the $`n=1`$ case was previously found by Rice (1954), and $`n=2`$ by Berry (1978)). For general $`n,`$ the density $`d_n`$ (1.2) is expressed as a probability integral with density function (2.1) and correlations given by (3.5),(3.6). Therefore

$`d_n=\delta ^n(\mathrm{v}_\mathrm{i}\mathbf{)}\mathbf{|}𝐝𝐞𝐭\mathrm{v}_{\mathrm{i}\mathbf{,}\mathrm{j}}\mathbf{|}\mathbf{}`$
$`={\displaystyle \frac{1}{(2\pi )^{n(n+1)/2}F_0^{n^2/2}}}{\displaystyle ^nv_i^{n^2}v_{i,j}\delta ^n(\mathrm{v}\mathbf{)}\mathbf{|}𝐝𝐞𝐭\mathrm{v}_{\mathrm{i}\mathbf{,}\mathrm{j}}\mathbf{|}\mathrm{𝐞𝐱𝐩}\mathbf{\left(}\mathbf{}\frac{1}{2}\underset{\mathrm{i}\mathbf{=}1}{\overset{\mathrm{n}}{\mathbf{}}}\mathrm{v}_\mathrm{i}^2\mathbf{}\frac{1}{2}\underset{\mathrm{i}\mathbf{,}\mathrm{j}\mathbf{=}1}{\overset{\mathrm{n}}{\mathbf{}}}\mathrm{v}_{\mathrm{i}\mathbf{,}\mathrm{j}}^2\mathbf{\right)}}`$
$`={\displaystyle \frac{F_0^{n/2}}{(2\pi )^{n(n+1)/2}}}{\displaystyle ^{n^2}v_{i,j}|detv_{i,j}|\mathrm{exp}\left(\frac{1}{2}\underset{i,j=1}{\overset{n}{}}v_{i,j}^2\right)},`$ (4.1)
where in the final line the $`\delta `$-functions in the $`v_i`$ were integrated, and the $`v_{i,j}`$ were rescaled (each by $`\sqrt{F_0}`$) to be dimensionless. The remaining integral is solved in B, giving the density of zeros in $`n`$ dimensions

$$d_n=(C_0^{\prime \prime })^{n/2}\frac{(\frac{n1}{2})!}{\pi ^{(n+1)/2}}=(C_0^{\prime \prime })^{n/2}\frac{(n1)!\sigma _{n1}}{(2\pi )^n}.$$
(4.2)
In this expression, $`\sigma _{n1}`$ is the surface area of the unit $`(n1)`$-sphere in $`n`$ dimensions, given by (2.4). As is common in such problems in statistical geometry, the result is a spectral quantity ($`(C_0^{\prime \prime })^{n/2}`$) multiplied by a geometric factor.

The scheme of the previous section may now be applied to calculate the topological charge correlation function for zeros in these gaussian random vector fields. Implementation of the scheme is facilitated by the fact that the components $`v_i`$ are completely independent, and $`m=n^2.`$ The submatrix $`𝐊_n`$ of the correlation matrix $`𝚺_n`$ only depends on the correlations of the components of the vectors $`\mathrm{v}_\mathrm{A}\mathbf{,}\mathrm{v}_\mathrm{B}\mathbf{;}`$ from section 3, the only correlations that do not vanish are $`v_{Ai}^2=v_{Bi}^2=1`$ and $`v_{Ai}v_{Bi}=C`$ for $`i=1,\mathrm{},n.`$ It is easy to see that

$$det𝐊_n=(1C^2)^n.$$
(4.3)
From (2.12), the other necessary ingredient of the correlation function scheme is the components of the matrix $`𝚵,`$ defined in (2.5). The elements of $`𝚵`$ are labelled by the multiindex of the components $`v_{Ai,j},v_{Bk,l};`$ using the correlations (3.5), (3.6) and the arguments of A, particularly equation (1.5), the only nonvanishing elements are

$`\mathrm{\Xi }_{(Ai,1)(Ai,1)}=\mathrm{\Xi }_{(Bi,1)(Bi,1)}=F_0E^2/(1C^2)`$
$`\mathrm{\Xi }_{(Ai,1)(Bi,1)}=\mathrm{\Xi }_{(Bi,1)(Ai,1)}=F_0E^2/(1C^2)`$
$`\mathrm{\Xi }_{(Ai,j)(Ai,j)}=\mathrm{\Xi }_{(Bi,j)(Bi,j)}=F_0`$
$`\mathrm{\Xi }_{(Ai,j)(Bi,j)}=\mathrm{\Xi }_{(Bi,j)(Ai,j)}=H\mathrm{for}i=1,\mathrm{},n,j=2,\mathrm{},n.`$ (4.4)
The problem remains to use these components and (2.12) to evaluate (2.11). By (1.1), (2.8),

$$𝒥_A=det_{Ai,j},$$
(4.5)
(similarly for $`𝒥_B`$). Each of the $`n!`$ summands in this determinant is an $`n`$-fold product $`\mathrm{sign}\sigma _i_{Ai,\sigma (i)},`$ where $`\sigma `$ is a permutation of $`1,\mathrm{},n.`$ Therefore

$$𝒥_A𝒥_B=\underset{\sigma ,\sigma ^{}}{\overset{\mathrm{permutations}}{}}\mathrm{sign}\sigma \mathrm{sign}\sigma ^{}(1)^n\underset{i,j=1}{\overset{n}{}}_{Ai,\sigma (i)}_{Bj,\sigma (i)}.$$
(4.6)
From (2.12), the result of one of these summands acting on $`\mathrm{exp}(\mathrm{t}\mathbf{}\mathrm{\Xi }\mathrm{t}\mathbf{/}2\mathbf{)}`$ and setting $`\mathrm{t}\mathbf{=}\mathrm{𝟎}`$ is nonzero if there is a pairing of these multiindices where the corresponding elements of $`𝚵`$ are nonzero. From (4.4), this is only the case when the permutations $`\sigma ,\sigma ^{}`$ are the same. Thus, from (2.12) and (4.4),

$`D_n`$ $`=`$ $`{\displaystyle \underset{\sigma }{\overset{\mathrm{permutations}}{}}}{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{\Xi }_{(Ai,\sigma (i))(Bi,\sigma (i))}`$ (4.7)
$`=`$ $`n!{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{\Xi }_{(Ai,i)(Bi,i)}`$
$`=`$ $`n!(FCE^2/(1C^2))H^{n1}.`$
This, together with (4.3), can now be put into (2.10) to give

$`g_n(r)`$ $`=`$ $`{\displaystyle \frac{n!}{2\pi d_n^2}}{\displaystyle \frac{(F(1C^2)CE^2)}{(1C^2)^{3/2}}}\left({\displaystyle \frac{H}{2\pi \sqrt{1C^2}}}\right)^{n1}`$ (4.8)
$`=`$ $`{\displaystyle \frac{n!}{2\pi d_n^2}}{\displaystyle \frac{(C^{\prime \prime }(1C^2)+CC^2)}{(1C^2)^{3/2}}}\left({\displaystyle \frac{C^{}}{2\pi r\sqrt{1C^2}}}\right)^{n1}`$
$`=`$ $`{\displaystyle \frac{(n1)!}{(2\pi )^nd_n^2r^{n1}}}{\displaystyle \frac{h_n(r)}{r}},`$
where the function $`h_n(r)`$ is defined

$$h_n(r)=\left(\frac{C^{}}{\sqrt{1C^2}}\right)^n.$$
(4.9)
This decays to 0 as $`r\mathrm{},`$ and when $`r=0,`$

$$h_n(0)=(C_0^{\prime \prime })^{n/2}$$
(4.10)
(the quantity $`C_0^{\prime \prime }=C^{\prime \prime }(0)`$ is always negative, since $`C`$ is the Fourier transform of the positive power spectrum $`\mathrm{\Pi }`$).

The topological charge correlation function of zeros of $`n`$-dimensional gaussian random vector fields $`g(r)`$ was first calculated by Halperin (1981), in a form equivalent to (4.8), but without proof. This function was also derived by different means by Liu and Mazenko (1992), and in the $`n=1`$ case by Rice (1954), and the $`n=2`$ case by Berry and Dennis (2000) and Foltin (2003a). $`g(r)`$ is plotted in figure 1 for two choices of the field correlation function $`C`$ for $`n=2.`$ When $`C(r)=J_0(r),`$ $`g_2`$ is oscillatory; when $`C(r)=\mathrm{exp}(r^2/2),`$ it increases monotonically to zero as $`r\mathrm{}.`$

## 5 Critical points in two dimensions

In this section, the scheme of section 2 is used to calculate the topological charge correlation function of critical points of isotropic gaussian random functions in the plane, that is the Poincaré index correlation function of isotropic random surfaces.

The gaussian random function examined shall be written $`f=f(\mathrm{r}\mathbf{)}\mathbf{=}\mathrm{f}\mathbf{(}\mathrm{x}\mathbf{,}\mathrm{y}\mathbf{)}`$ (where $`x=r_1,y=r_2`$), which will be assumed stationary and isotropic, so the expressions in section 3 may be used. In particular, $`f`$ and its first derivatives have the correlations (3.5), (3.6), as well as further correlations involving second derivatives, described below. As in the previous section, for convenience in calculations, the two points $`A`$ and $`B`$ are separated only in the $`x`$ coordinate.

At a critical point, the gradient $`f=(f_{,x},f_{,y})`$ is zero. The critical point jacobian $`𝒥_\mathrm{c}`$ whose sign defines the topological charge (Poincaré index) is the hessian determinant

$$𝒥_\mathrm{c}=detf_{,ij}=f_{,xx}f_{,yy}f_{,xy}^2;$$
(5.1)
this is the gaussian curvature of the surface. Unlike a general 2-dimensional random vector field (as in the previous section), the gradient field $`f`$ is irrotational, which gives rise to relationships and correlations between the derivatives of the components (e.g. $`f_{,xy}=f_{,yx}`$, whereas before, $`v_{1,2}v_{2,1}`$ in general). This makes the explicit computation of the topological charge correlation function more difficult than in the previous section, even in the 2-dimensional case that is considered here; the scheme of section 2 applies to gradient zeros in fields of any dimension.

The statistical properties of critical points of gaussian random fields in 2 dimensions were considered by Longuet-Higgins (1957a,b); he found that the density $`d_\mathrm{c}`$ of critical points is (Longuet-Higgins 1957b equations (71), (78)):

$$d_\mathrm{c}=\left|\frac{2C_0^{(4)}}{3\pi \sqrt{3}C_0^{\prime \prime }}\right|,$$
(5.2)
($`C_0^{(4)}`$ denotes the fourth derivative of $`C,`$ evaluated at $`r=0`$). The density of saddles equals the density of extrema (maxima and minima), and the density of maxima equals the density of minima. The probability density function of the gaussian curvature $`𝒥`$, despite its asymmetry (Longuet-Higgins 1958 equation (7.14), Dennis 2002 equation (57)), has zero first moment, implying the average topological charge $`\delta ^2(f)𝒥_\mathrm{c}`$ is zero, as expected.

The topological charge correlation function is again calculated using the scheme of section 2, particularly (2.10). Therefore, the vector $`\mathrm{u}`$ (2.3) of gaussian random variables, in a convenient ordering, is

$$\mathrm{u}_𝐜\mathbf{=}\mathbf{(}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{y}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{y}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{y}}\mathbf{)}$$
(5.3)
with correlation matrix (c.f. (2.2))

$$𝚺_\mathrm{c}=\left(\begin{array}{cccccccccc}M_0& L_0& M& L& 0& 0& 0& G& 0& 0\\ L_0& M_0& L& N& 0& 0& 0& I& 0& 0\\ M& L& M_0& L_0& 0& 0& G& 0& 0& 0\\ L& N& L_0& M_0& 0& 0& I& 0& 0& 0\\ 0& 0& 0& 0& L_0& L& 0& 0& 0& I\\ 0& 0& 0& 0& L& L_0& 0& 0& I& 0\\ 0& 0& G& I& 0& 0& F_0& F& 0& 0\\ G& I& 0& 0& 0& 0& F& F_0& 0& 0\\ 0& 0& 0& 0& 0& I& 0& 0& F_0& H\\ 0& 0& 0& 0& I& 0& 0& 0& H& F_0\end{array}\right),$$
(5.4)
where the correlations between elements of $`\mathrm{u}_𝐜`$ are computed to be

$`G`$ $`\left[_x^3C\right]_{x=r,y=0}`$ $`=C^{(3)}`$
$`I`$ $`\left[_x^2_yC\right]_{x=r,y=0}`$ $`=(rC^{\prime \prime }C^{})/r^2`$
$`L`$ $`\left[_x^2_y^2C\right]_{x=r,y=0}`$ $`=(r^2C^{(3)}2rC^{\prime \prime }+2C^{})/r^3`$
$`M`$ $`\left[_x^4C\right]_{x=r,y=0}`$ $`=C^{(4)}`$
$`N`$ $`\left[_y^4C\right]_{x=r,y=0}`$ $`=3(rC^{\prime \prime }C^{})/r^3`$
$`G_0`$ $`=I_0=0,`$ $`M_0=N_0=3L_0=C_0^{(4)}.`$ (5.5)
The last line gives the special value of these correlations when $`r=0.`$

The matrix $`𝐊_\mathrm{c}`$ is the $`4\times 4`$ lower right submatrix of $`𝚺_\mathrm{c},`$ and has determinant

$$det𝐊_\mathrm{c}=(F_0^2H^2)(F_0^2F^2).$$
(5.6)
The pair of differential jacobian operators are, from (5.1),

$$𝒥_A𝒥_B=_1_2_3_4+_5^2_6^2_1_2_6^2_3_4_5^2.$$
(5.7)
Using (2.12), the result of these operators acting is on $`\mathrm{exp}(\mathrm{t}\mathbf{}\mathrm{\Xi }\mathbf{}\mathrm{t}\mathbf{/}2\mathbf{)}`$ and setting $`\mathrm{t}\mathbf{=}\mathrm{𝟎}`$ (c.f. (2.10),(2.11)) is

$$D_\mathrm{c}=\mathrm{\Xi }_{14}\mathrm{\Xi }_{23}+\mathrm{\Xi }_{13}\mathrm{\Xi }_{24}2\mathrm{\Xi }_{16}\mathrm{\Xi }_{26}+\mathrm{\Xi }_{12}\mathrm{\Xi }_{34}2\mathrm{\Xi }_{35}\mathrm{\Xi }_{45}\mathrm{\Xi }_{34}\mathrm{\Xi }_{55}+2\mathrm{\Xi }_{56}^2\mathrm{\Xi }_{12}\mathrm{\Xi }_{66}+\mathrm{\Xi }_{55}\mathrm{\Xi }_{66},$$
(5.8)
where the necessary entries of $`𝚵_\mathrm{c},`$ are found using Jacobi’s determinant theorem in (1.5); as an example,

$$\mathrm{\Xi }_{24}=M_0FI^2/(F_0^2F^2).$$
(5.9)
The topological charge correlation function $`g_\mathrm{c}(r)`$ for critical points is obtained by substituting (5.6), (5.8) (with all terms like (5.9) found using (1.5) into (2.10)). This expression is complicated and not very illuminating, and is not given here. Upon substituting (3.5), (3.6), (5.5) in, one finds that $`g_\mathrm{c}`$ can be written as a perfect derivative (c.f. (4.8)),

$$g_\mathrm{c}(r)=\frac{1}{4\pi ^2d_\mathrm{c}^2r}\frac{h_\mathrm{c}}{r},$$
(5.10)
where

$`h_\mathrm{c}(r)={\displaystyle \frac{(C^{\prime \prime }C^{}/r)}{r\sqrt{(C_0^{\prime \prime 2}C^{\prime \prime 2})(C_0^{\prime \prime 2}C^2/r^2)}}}[{\displaystyle \frac{C^{(3)}(3C_0^{\prime \prime 2}2C^{\prime \prime 2}C^{\prime \prime }C^{}/r)}{C_0^{\prime \prime 2}C^{\prime \prime 2}}}`$
$`\text{ }{\displaystyle \frac{(C^{\prime \prime }C^{}/r)(3C_0^{\prime \prime 2}2C^2/r^2C^{\prime \prime }C^{}/r)}{r(C_0^{\prime \prime 2}C^2/r^2)}}]`$
$`={\displaystyle \frac{1}{r(C^{\prime \prime }C^{}/r)}}{\displaystyle \frac{}{r}}{\displaystyle \frac{(C^{\prime \prime }C^{}/r)^3}{\sqrt{(C_0^{\prime \prime 2}C^{\prime \prime 2})(C_0^{\prime \prime 2}C^2/r^2)}}}.`$ (5.11)
This process of finding $`h_\mathrm{c}`$ is long and tedious, and details are omitted here. It is easy to see that $`h_\mathrm{c}0`$ as $`r\mathrm{};`$ it is straightforward, by Taylor expanding derivatives of $`C,`$ to show that

$$h_\mathrm{c}(0)=\frac{4C_0^{(4)}}{3\sqrt{3}C_0^{\prime \prime }}=2\pi d_\mathrm{c}.$$
(5.12)
The critical point topological charge correlation function for two particular field correlation functions is shown in figure 2. As with the 2-dimensional vector case plotted in figure 1, the properties of the correlation function (oscillatory, exponential decay, etc) are similar to that of the underlying field correlation function $`C(r),`$ on which the correlation function depends. Features of interest in these plots include the sharp initial minimum of (a), and the fact that the monotonic decay in (b) is from above, not below (in contrast to its counterpart in figure 1. Nevertheless, the form of $`h_\mathrm{c}`$ is significantly more complicated than $`h_n,`$ especially when $`n=2.`$

## 6 Umbilic points

Less well-known than critical points, umbilic points are geometric point singularity features associated with the second derivative of $`f`$ \- namely where hessian matrix of second derivatives $`_{ij}f`$ is degenerate (Berry and Hannay 1977, Porteous 2001, Hilbert and Cohn-Vossen 1952). Geometrically, the principal axes of gaussian curvature are not defined at these points. The eigenvalues of the hessian coincide when $`f_{,xx}=f_{,yy},f_{,xy}=0.`$ Umbilic points are therefore zeros of the 2-dimensional vector field

$$\mathrm{v}_𝐮\mathbf{=}\mathbf{(}\mathbf{(}\mathrm{f}_{\mathbf{,}\mathrm{x}\mathrm{x}}\mathbf{}\mathrm{f}_{\mathbf{,}\mathrm{y}\mathrm{y}}\mathbf{)}\mathbf{/}2\mathbf{,}\mathrm{f}_{\mathbf{,}\mathrm{x}\mathrm{y}}\mathbf{)}\mathbf{.}$$
(6.1)
The factor of half in the first term ensures that $`\mathrm{v}_𝐮`$ is statistically rotation invariant.

An umbilic point has an index, determined geometrically by the sense of rotation of the principal axes of curvature around the umbilic point, and the index is generically $`\pm 1/2`$ (Berry and Hannay 1977); only the sign of the index is of interest here, and this is determined by the appropriate jacobian $`𝒥_\mathrm{u}`$ on $`\mathrm{v}_𝐮`$,

$$2𝒥_\mathrm{u}=f_{,xxx}f_{,xyy}+f_{,yyy}f_{,xxy}f_{,xyy}^2f_{,xxy}^2,$$
(6.2)
depending on the third partial derivatives of $`f.`$ The calculation of the topological charge correlation function for umbilic points can proceed according to the scheme of section 2, in a similar way to the corresponding calculation for critical points.

Umbilic points for isotropic random functions was considered by Berry and Hannay; they found the density $`d_\mathrm{u}`$ to be (Berry and Hannay 1977 equation (34))

$$d_\mathrm{u}=\left|\frac{3C_0^{(6)}}{10\pi C_0^{(4)}}\right|,$$
(6.3)
($`C_0^{(6)}`$ is the sixth derivative of $`C`$ at 0) and that the average index $`q`$ is zero (the separate densities are $`0.5d_\mathrm{u}`$ for stars, $`0.05279d_\mathrm{u}`$ for monstars, and $`0.44721d_\mathrm{u}`$ for lemons). In the present work, the distinction between monstars and lemons, which both have positive index, is not used.

The ordering of the vector of gaussian random functions $`\mathrm{u}_𝐮`$ is chosen

$`\mathrm{u}_𝐮\mathbf{=}\mathbf{(}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}\mathrm{x}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}\mathrm{y}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}\mathrm{x}\mathrm{x}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}\mathrm{y}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{x}\mathrm{x}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{A}\mathbf{,}\mathrm{y}\mathrm{y}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{x}\mathrm{x}\mathrm{y}}\mathbf{,}\mathrm{f}_{\mathrm{B}\mathbf{,}\mathrm{y}\mathrm{y}\mathrm{y}}\mathbf{,}`$
$`(f_{A,xx}f_{A,yy})/2,(f_{B,xx}f_{B,yy})/2,f_{A,xy},f_{B,xy}).`$ (6.4)
The correlation matrix (2.2) is

$$𝚺_\mathrm{u}=\left(\begin{array}{cccccccccccc}S_0& T_0& S& T& 0& 0& 0& 0& 0& X& 0& 0\\ T_0& T_0& T& U& 0& 0& 0& 0& 0& Y& 0& 0\\ S& T& S_0& T_0& 0& 0& 0& 0& X& 0& 0& 0\\ T& U& T_0& T_0& 0& 0& 0& 0& Y& 0& 0& 0\\ 0& 0& 0& 0& T_0& T_0& T& U& 0& 0& 0& Q\\ 0& 0& 0& 0& T_0& S_0& U& V& 0& 0& 0& R\\ 0& 0& 0& 0& T& U& T_0& T_0& 0& 0& Q& 0\\ 0& 0& 0& 0& U& V& T_0& S_0& 0& 0& R& 0\\ 0& 0& X& Y& 0& 0& 0& 0& L_0& W& 0& 0\\ X& Y& 0& 0& 0& 0& 0& 0& W& L_0& 0& 0\\ 0& 0& 0& 0& 0& 0& Q& R& 0& 0& L_0& L\\ 0& 0& 0& 0& Q& R& 0& 0& 0& 0& L& L_0\end{array}\right),$$
(6.5)
where $`W(M+N2L)/4,X(PQ)/2,Y(QR)/2,`$ and

$`P`$ $`\left[_x^5C\right]_{x=r,y=0}`$ $`=C^{(5)},`$
$`Q`$ $`\left[_x^3_y^2C\right]_{x=r,y=0}`$ $`=(r^3C^{(4)}3r^2C^{(3)}+6rC^{\prime \prime }6C^{})/r^4,`$
$`R`$ $`\left[_x_y^4C\right]_{x=r,y=0}`$ $`=3(r^2C^{(3)}3rC^{\prime \prime }+3C^{})/r^4,`$
$`S`$ $`\left[_x^6C\right]_{x=r,y=0}`$ $`=C^{(6)},`$
$`T`$ $`\left[_x^4_y^2C\right]_{x=r,y=0}`$ $`=(r^4C^{(5)}4r^3C^{(4)}+12r^2C^{(3)}24rC^{\prime \prime }+24C^{})/r^5,`$
$`U`$ $`\left[_x^2_y^4C\right]_{x=r,y=0}`$ $`=3(r^3C^{(4)}5r^2C^{(3)}+12rC^{\prime \prime }12C^{})/r^5,`$
$`V`$ $`\left[_y^6C\right]_{x=r,y=0}`$ $`=15(r^2C^{(3)}3rC^{\prime \prime }+3C^{})/r^5,`$
$`P_0`$ $`=Q_0=R_0=0,`$ $`S_0=5T_0=5U_0=V_0=C_0^{(6)}.`$ (6.6)
The matrix $`𝐊_\mathrm{u},`$ defined using (2.5), has determinant

$$det𝐊_\mathrm{u}=(L_0^2(M+N2L)^2/16)(L_0^2L^2).$$
(6.7)
The result of the jacobian derivative operators (2.11) gives

$`D_\mathrm{u}=\mathrm{\Xi }_{12}^2+\mathrm{\Xi }_{14}^2+\mathrm{\Xi }_{22}^2+2\mathrm{\Xi }_{24}^22\mathrm{\Xi }_{12}\mathrm{\Xi }_{22}+\mathrm{\Xi }_{13}\mathrm{\Xi }_{24}4\mathrm{\Xi }_{14}\mathrm{\Xi }_{24}2\mathrm{\Xi }_{12}\mathrm{\Xi }_{55}+2\mathrm{\Xi }_{12}\mathrm{\Xi }_{56}`$
$`2\mathrm{\Xi }_{22}\mathrm{\Xi }_{56}+\mathrm{\Xi }_{55}^2+\mathrm{\Xi }_{56}^2+2\mathrm{\Xi }_{57}^2+\mathrm{\Xi }_{58}^22\mathrm{\Xi }_{55}\mathrm{\Xi }_{56}4\mathrm{\Xi }_{57}\mathrm{\Xi }_{58}+\mathrm{\Xi }_{57}\mathrm{\Xi }_{68}.`$ (6.8)
The necessary entries of $`𝚵_\mathrm{u}`$ are found using (1.5). The resulting expression for $`D_\mathrm{u},`$ and therefore for $`g_\mathrm{u},`$ is very complicated, but may be reduced to the following form:

$$g_\mathrm{u}(r)=\frac{1}{4\pi ^2d_\mathrm{u}^2r}\frac{h_\mathrm{u}}{r},$$
(6.9)
where

$$h_\mathrm{u}(r)=\frac{}{r}\left(\frac{r(QR)^2}{4\sqrt{det𝐊_\mathrm{u}}}\right)+\frac{Q(P+R2Q)}{4\sqrt{det𝐊_\mathrm{u}}}.$$
(6.10)
It can be shown that $`h_\mathrm{u}(0)=2\pi d_\mathrm{u},`$ and $`h0`$ as $`r\mathrm{}.`$ $`g_\mathrm{u}(r)`$ is plotted in figure 3 for $`C(r)=J_0(r)`$ and $`\mathrm{exp}(r^2/2).`$ The most obvious feature of these two plots, compared to figures 1, 2, is that they have a negative maximum, a property that seems to be general for $`g_\mathrm{u},`$ although this has not been proved. It is unclear what the physical significance of this kink should be; mathematically, it probably arises from interference between the two summands in $`h_\mathrm{u}(r)`$ in (6.10).

## 7 Topological charge screening

Three particular charge correlation functions have been derived exactly (equations (4.8), (4.9), (5.10), (5.11), (6.9),(6.10)). In each case, the charge correlation function has the form

$$g(r)=\frac{(n1)!}{(2\pi )^nd^2r^{n1}}\frac{h(r)}{r},$$
(7.1)
where $`n`$ is the number of dimensions, $`d`$ is the density of zeros and $`h(r)`$ is a function such that

$$h(0)=\frac{(2\pi )^nd}{\sigma _{n1}(n1)!},h\mathrm{}\mathrm{as}r\mathrm{}.$$
(7.2)
The total charge density around a given (say positive) topological charge is therefore

$`d{\displaystyle ^n\mathrm{r}\mathrm{g}\mathbf{(}\mathrm{r}\mathbf{)}}`$ $`=`$ $`{\displaystyle \frac{(n1)!}{(2\pi )^nd}}\sigma _{n1}{\displaystyle _0^{\mathrm{}}}r{\displaystyle \frac{h(r)}{r}}`$ (7.3)
$`=`$ $`{\displaystyle \frac{(n1)!\sigma _{n1}}{(2\pi )^nd^2}}(h(\mathrm{})h(0))`$
$`=`$ $`1,`$
where the hypersphere area $`\sigma _{n1}`$ appears in the first line from integration in polar coordinates. It implies that the distribution of topological charges is such that every topological charge tends to be surrounded by charges of the opposite sign, such that the topological charge is ‘screened’ at large distances. This fact was noticed in the random vector case by Halperin (1981) and Liu and Mazenko (1992) (although the general zero density $`d_n`$ had not been determined explicitly) and is independent of the field correlation function $`C.`$ The derivation here shows that this is a more general phenomenon, possibly a universal feature of topological charge correlations for gaussian random fields. It should be noted that (7.3) is not necessarily satisfied for an arbitrary distribution of signed points; for instance, $`g(r)=0`$ always for Poisson points, for which there is no screening. $`\delta `$-function correlations at the origin are ignored in the following.

An analogy may be drawn from the theory of ionic liquids (Hansen and McDonald 1986); in a fluid or plasma, consisting of two species identical except for opposite (Coulomb) charges, the following Stillinger-Lovett sum rules (Stillinger and Lovett 1968a,b) are found to hold:

$`d{\displaystyle ^n\mathrm{r}\mathrm{g}\mathbf{(}\mathrm{r}\mathbf{)}}`$ $`=`$ $`1,`$ (7.4)
$`d{\displaystyle ^n\mathrm{r}\mathrm{r}^2\mathrm{g}\mathbf{(}\mathrm{r}\mathbf{)}}`$ $`=`$ $`a_n\mathrm{\Lambda }^2.`$ (7.5)
Here, $`g`$ is the charge-charge correlation function, $`\mathrm{\Lambda }`$ is a characteristic screening length dependent on temperature and $`n,`$ and $`a_n`$ is a constant dependent on dimensionality. These rules are discussed for $`n=3`$ by Hansen and McDonald (1986), Stillinger and Lovett 1968a,b, and $`n=2`$ by Jancovici 1987, Jancovici et al1994.

The screening relation (7.3) is equivalent to the first Stillinger-Lovett sum rule (7.4), which is derived using the statistical mechanics of pairwise, Coulomb interacting fluids. It is unclear whether the fact that topological and coulombic charges screen in the same way is coincidence, or evidence of some deeper connection between the two statistical theories.

It is natural to ask whether the topological charge correlation functions satisfy the second sum rule, which (upon integrating the left hand side of (7.5) by parts), depends on the integral of $`rh.`$ For $`n=2,`$ it was found (Berry and Dennis 2000) that for certain choices of $`C,`$ this integral may diverge. The slowest decay comes when $`C=J_0,`$ and by (4.9),

$$h_2(r)_r\mathrm{}\mathrm{cos}(r+\pi /4)^2/r\mathrm{when}C=J_0,$$
(7.6)
giving a logarithmic divergence for the second moment. For comparison, the critical and umbilic functions $`h_\mathrm{c},h_\mathrm{u},`$ (equations (5.11), (6.10) respectively) both give, for the same choice of $`C,`$

$$h_\mathrm{c}(r),h_\mathrm{u}(r)_r\mathrm{}\mathrm{cos}(2r)/r^2\mathrm{when}C=J_0,$$
(7.7)
implying that the second sum rule is satisfied generally for critical and umbilic points, although, as for the cases of random vectors where the integral converges, the screening length $`\mathrm{\Lambda },`$ defined in analogy to (7.5), depends on the choice of $`C.`$ For random vectors with $`n>2,`$ the second moment of $`g`$ always converges, because of the higher power of $`C^{}`$ appearing in (4.9) (also, the decay of $`C`$ will be faster, as discussed at the end of section 3).

Comparison may be drawn to the electrostatic analogy in random matrix theory, particularly in the case of the so-called Ginibre ensemble of $`N\times N`$ matrices whose entries are independent, identically distributed circular gaussian random variables (Ginibre 1965). The eigenvalues of these matrices are found to have exactly the same statistical behaviour as a 1-component 2-dimensional Coulomb gas of $`N`$ charges in a harmonic oscillator potential, and the 2-point density correlation function screens against a uniform background (i.e. $`d^2z(g_{\mathrm{Gin}}1)=1`$), and have finite second moment. Certain random polynomial analogues have zeros that can also be expressed as 2-dimensional Coulomb gases with additional interactions (Hannay 1998, Forrester and Honner 1999). The eigenvalues of random matrices (which may be expressed as the zeros of the characteristic polynomial) and zeros of random polynomials are all of the same sign, since they are zeros of complex analytic functions, and the density correlation functions of zeros are unique (there is no analogue of $`C`$).

There may be a danger in taking the analogy with fluids too far; for instance, the oscillations of the functions in figures 1,2,3a are reminiscent of those of charged fluids (e.g. Hansen and McDonald 1986). However, the physical causes for these oscillations are very different. In fluids, the oscillations usually arise from packing considerations (the ions themselves are of finite size, fixing the lengthscale, although in plasmas they are usually represented as point charges (Baus and Hansen 1980)). Topological charges, on the other hand, are points, and the oscillations in these figures originate from the oscillations in the underlying field correlation function $`C(r),`$ which is $`J_0`$ in this case. Although Halperin (1981) discusses the similarity between the short-range behaviour of the 2-dimensional vector correlation function (4.8), (4.9) and Kosterlitz-Thouless theory, the present situation is more general, both in that the (possibly long-range) screening is exact, and that the results hold for any reasonable field correlation function $`C(r).`$

## 8 Discussion

Using a general scheme for calculating topological charge correlation functions, three particular correlation functions were found explicitly, and were found always to satisfy a screening relation (7.3).

The scheme of section 2 used to calculate the charge correlation functions is very general, and can be generalised to calculate the charge-charge correlation function between two different sets of topological charges - for instance with a critical point at $`A,`$ and an umbilic point at $`B.`$ Although not done so here, the scheme may be applied to anisotropic fields.

It is rather more difficult to calculate the density correlation function for topological charges (the analogue of (1.3) where the moduli of the jacobians are taken). It was calculated for 1-dimensional vectors (i.e. random functions in 1 dimension) by Rice (1954), and for 2-dimensional vectors (realised as complex scalars) by Berry and Dennis (2000), Saichev et al(2001). However, it has not been possible to generalise such methods to the case of critical points. Also, numerical evidence (Freund and Wilkinson 1998) suggests that the correlation function of extrema signed by the sign of their laplacian ($`+1`$ for minima, $`1`$ for maxima) also satisfies the screening relation (7.4). All such functions would be needed to calculate the partial correlation functions between the species (e.g. maxima with maxima, maxima with saddles, maxima with minima, etc), which would give a complete statistical picture.

The presence of boundaries in the random function will affect the statistical properties of topological charges (e.g. Berry (2002) for nodal points in the plane), and it is possible that there may be some further analogy with the physics of interfaces of Coulomb fluids. The scheme of section 2 ought to be adaptable to calculate charge correlation functions in this case.

Only zeros of fields linear in gaussian random functions have been considered here, although the density of others may be calculated, for instance, in addition to nodal points, a 2-dimensional gaussian random complex scalar has critical points of its modulus squared (Weinrib and Halperin 1982) and its argument (Dennis 2001a). The scheme employed here cannot be used to calculate correlation functions for these, although numerical evidence (I Freund, personal communication) suggests that the critical points of argument (together with the nodal points) do screen, therefore adding to the cases shown here. It is tempting to conjecture that topological charge screening may be a universal phenomenon in gaussian random fields.

I would like to thank Michael Berry and Robert Evans for useful discussions, John Hannay for discussions leading to the arguments in Appendix B, and Isaac Freund for correspondence. This research was supported by the Leverhulme Trust.

## Appendix A Jacobi’s determinant theorem

Let $`𝐀`$ be a square matrix. The minor $`_{j_1\mathrm{}j_k}^{i_1\mathrm{}i_k}(𝐀)`$ is the determinant of the $`k\times k`$ submatrix of $`𝐀`$ with rows $`i_1,\mathrm{},i_k`$ and columns $`j_1,\mathrm{},j_k.`$ $`\overline{}_{j_1,\mathrm{},j_k}^{i_1,\mathrm{},i_k}(𝐀)`$ shall be used to denote the complementary minor, that is, the minor of the submatrix of $`𝐀`$ with rows $`i_1,\mathrm{},i_k`$ and columns $`j_1,\mathrm{},j_n`$ excluded. Then Jacobi’s determinant theorem (Jeffreys and Jeffreys 1956, page 135) states

$$det𝐀_{j_1,\mathrm{},j_k}^{i_1,\mathrm{},i_k}(𝐀^1)=(1)^{i_1+\mathrm{}+i_k+j_1+\mathrm{}+j_k}\overline{}_{j_1,\mathrm{},j_k}^{i_1,\mathrm{},i_k}(𝐀).$$
(1.1)
Applying this to $`𝚺`$ in (2.5), and choosing $`𝚵^1`$ as the submatrix whose determinant is the minor of $`𝚺^1,`$

$`det𝚺det𝚵^1`$ $`=`$ $`det𝚺_{1,\mathrm{},2m}^{1,\mathrm{},2m}(𝚺^1)`$ (1.2)
$`=`$ $`(1)^{1+\mathrm{}+2m+1+\mathrm{}+2m}_{2m+1,\mathrm{},2(m+n)}^{2m+1,\mathrm{},2(m+n)}(𝚺)`$
$`=`$ $`det𝐊,`$
from which (2.6) follows directly.

Jacobi’s theorem can also be used to find the elements $`\mathrm{\Xi }_{ij}`$ of the inverse reduced inverse matrix $`𝚵`$ in (2.5), needed for (2.12). In this case, (1.1) is applied twice, once on the matrix $`𝚵^1,`$ and once on $`𝚺.`$ Therefore

$`\mathrm{\Xi }_{ij}`$ $`=`$ $`_j^i(𝚵)`$ (1.5)
$`=`$ $`(1)^{i+j}\overline{}_j^i(𝚵^1)/det𝚵^1`$
$`=`$ $`(1)^{i+j}det𝚵_{1,\mathrm{},j1,j+1,\mathrm{},2m}^{1,\mathrm{},i1,i+1,\mathrm{},2m}(𝚺^1)`$
$`=`$ $`det𝚵_{j,2m+1,\mathrm{},2m+2n}^{i,2m+1,\mathrm{},2m+2n}(𝚺)`$
$`=`$ $`_{j,2m+1,\mathrm{},2m+2n}^{i,2m+1,\mathrm{},2m+2n)}(𝚺)/det𝐊,`$
$`=`$ $`det𝐊^1det\left(\begin{array}{cc}\mathrm{\Sigma }_{ij}& \\ & 𝐊\end{array}\right)`$
where in the last line $``$ represents the terms $`\mathrm{\Sigma }_{i,k},\mathrm{\Sigma }_{k,j}`$ where $`k=2m+1,\mathrm{},2m+2n.`$ Thus the $`\mathrm{\Xi }_{ij}`$ appearing in the expression for $`D`$ in (2.12), is the determinant of the $`(2n+1)\times (2n+1)`$ submatrix comprised of the $`i`$th row and $`j`$th column of $`𝚺,`$ and the submatrix $`𝐊.`$

## Appendix B Calculation of the density of vector zeros (4.1)

In order to integrate (4.1), the following must be integrated

$$𝒱=^{n^2}v_{i,j}|detv_{i,j}|\mathrm{exp}\left(\frac{1}{2}\underset{i,j=1}{\overset{n}{}}v_{i,j}^2\right).$$
(2.1)
This is, mathematically, the average hypervolume of an $`n`$-dimensional parallelepiped specified by gaussian random vectors $`\mathrm{w}_1\mathbf{=}\mathbf{(}\mathrm{v}_{1\mathbf{,}1}\mathbf{,}\mathbf{}\mathrm{v}_{1\mathbf{,}\mathrm{n}}\mathbf{)}\mathbf{,}\mathbf{}\mathbf{,}\mathrm{w}_\mathrm{n}\mathbf{=}\mathbf{(}\mathrm{v}_{\mathrm{n}\mathbf{,}1}\mathbf{,}\mathbf{}\mathbf{,}\mathrm{v}_{\mathrm{n}\mathbf{,}\mathrm{n}}\mathbf{)}\mathbf{.}`$ These gaussian random vectors are identically distributed isotropically in $`n`$-dimensional space. The hypervolume is nonzero exactly when the set of $`n`$ vectors is linearly independent.

This hypervolume may be found explicitly in a manner reminiscent of the Gram-Schmidt orthogonalization procedure for vectors. Geometrically,

$`\mathrm{volume}\mathrm{of}\mathrm{parallelepiped}=\mathrm{length}\mathrm{of}\mathrm{w}_1\mathbf{\times }\mathrm{𝐥𝐞𝐧𝐠𝐭𝐡}\mathrm{𝐨𝐟}\mathrm{w}_2\mathrm{𝐨𝐫𝐭𝐡𝐨𝐠𝐨𝐧𝐚𝐥}\mathrm{𝐭𝐨}\mathrm{w}_1\mathbf{\times }\mathbf{}`$
$`\times \mathrm{length}\mathrm{of}\mathrm{w}_\mathrm{n}\mathrm{𝐨𝐫𝐭𝐡𝐨𝐠𝐨𝐧𝐚𝐥}\mathrm{𝐭𝐨}\mathrm{𝐬𝐩𝐚𝐧}\mathbf{\{}\mathrm{w}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{w}_{\mathrm{n}\mathbf{}1}\mathbf{\}}\mathbf{.}`$ (2.2)
A given factor in this product is therefore the average length of the gaussian random vector $`\mathrm{w}_\mathrm{j}`$ in the orthogonal complement of a $`(j1)`$-dimensional subspace $`\mathrm{span}\{\mathrm{w}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{w}_{\mathrm{j}\mathbf{}1}\mathbf{\}}\mathbf{.}`$

Now, since the vector $`\mathrm{w}_\mathrm{j}`$ is isotropic, it may be represented identically in any choice of orthonormal basis of $`n`$-dimensional space; in particular, its first $`j1`$ components $`v_{j,1},\mathrm{},v_{j,j1}`$ may be chosen to be in $`\mathrm{span}\{\mathrm{w}_1\mathbf{,}\mathbf{}\mathbf{,}\mathrm{w}_{\mathrm{j}\mathbf{}1}\mathbf{\}}`$ (as in the Gram-Schmidt procedure). The total contribution of $`\mathrm{w}_\mathrm{j}`$ to the integral in (4.1) involves the average length of the vector made up of the other components $`v_{j,j},\mathrm{},v_{j,n}.`$ Where $`k=nj1,`$ this is

$`{\displaystyle ^n\mathrm{w}_\mathrm{j}\sqrt{\mathrm{v}_{\mathrm{j}\mathbf{,}1}^2\mathbf{+}\mathbf{}\mathbf{+}\mathrm{v}_{\mathrm{j}\mathbf{,}\mathrm{n}}^2}\mathrm{𝐞𝐱𝐩}\mathbf{\left(}\mathbf{}\frac{1}{2}\underset{\mathrm{i}\mathbf{=}1}{\overset{\mathrm{n}}{\mathbf{}}}\mathrm{v}_{\mathrm{j}\mathbf{,}\mathrm{i}}^2\mathbf{\right)}}`$
$`=\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}v_{j,1}\mathrm{exp}(v_{j,1}^2/2)\right]^{j1}{\displaystyle _{𝒮_{k1}}}^{k1}\mathrm{\Omega }_{k1}{\displaystyle _0^{\mathrm{}}}\rho \rho ^k\mathrm{exp}(\rho ^2/2)`$ (2.3)
where, in the second line, the first $`m1`$ components have been pulled out as trivial gaussians, integrating to $`(2\pi )^{(j1)/2}=(2\pi )^{(nk)/2};`$ the remaining $`k`$ integrals are the average length of a gaussian random vector in $`k`$-dimensional space. This integral has been converted to polar coordinates, with $`\omega _{k1}`$ the solid angle infinitesimal on the unit $`(k1)`$-sphere $`𝒮_{k1},`$ and $`\rho `$ is the radius. It is well known that the surface area $`\sigma _{k1}`$ of the unit $`(k1)`$-hypersphere is

$$\sigma _{k1}=^{k1}\mathrm{\Omega }_{k1}=\frac{2\pi ^{k/2}}{(\frac{k2}{2})!}.$$
(2.4)
The $`\rho `$ integral in (2.3) is $`2^{(k1)/2}((k1)/2)!.`$ Therefore, the numerical part of (4.1) is $`1/(2\pi )^{n(n+1)/2}`$ times the product in (2.2), with each term in the product, now labelled by $`k,`$ given by the expression (2.3). Therefore

$`𝒱`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^{n(n+1)/2}}}{\displaystyle \underset{k=1}{\overset{n}{}}}(2\pi )^{(nk)/2}\times {\displaystyle \frac{2\pi ^{k/2}}{(\frac{k2}{2})!}}\times 2^{(k1)/2}\left({\displaystyle \frac{k1}{2}}\right)!`$ (2.5)
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^{n(n+1)/2}}}{\displaystyle \underset{k=1}{\overset{n}{}}}(2\pi )^{(n+1)/2}\pi ^{n/2}{\displaystyle \frac{(\frac{k1}{2})!}{(\frac{k2}{2})!}}`$
$`=`$ $`{\displaystyle \frac{(\frac{n1}{2})!}{\pi ^{(n+1)/2}}}.`$
This value agrees with that stated by Halperin (1981), Liu and Mazenko (1992) of $`1/\pi `$ ($`n=1`$), $`1/2\pi `$ ($`n=2`$) and $`1/\pi ^2`$ ($`n=3`$).

## References


format: text

text:
# Hidden Polynomial(s) Cryptosystems

## 1. Introduction

This paper focuses on Hidden Monomial Cryptosystems, a class of public key cryptosystems first proposed by Imai and Matsumoto \[imai1\]. In this class, the public key is a set of polynomial nonlinear equations. The private key is the set of parameters that the user chooses to construct the equations. Before we discuss our variation, we review briefly a simplified version of the original cryptosystem, better described in \[koblitz\]. The characters met throughout this paper are:

* Alice who wants to receive secure messages;
* Bob who wants to send her secure messages;
* Eve, the eavesdropper.

Alice takes two finite fields $`𝔽_q<𝕂`$, $`q`$ a power of $`2`$, and $`\beta _1,\beta _2,\mathrm{},\beta _n`$ a basis of $`𝕂`$ as an $`𝔽_q`$-vector space. Next she takes $`0<h<q^n`$ such that $`h=q^\theta +1`$, and $`gcd(h,q^n1)=1`$. Then she takes two generic vectors $`𝐮=(u_1,\mathrm{},u_n)`$ and $`𝐯=(v_1,\mathrm{},v_n)`$ upon $`𝔽_q`$, and puts111In this paper we reserve boldface to the elements of $`𝕂`$ thought as vectors upon $`𝔽_q`$ in the fixed private basis. They are considered vectors or field elements, as convenient, without further notice. This shift in practice takes a Chinese Remainder Theorem. In order to avoid boring repetitions, Cryptosystem and Scheme are used like synonyms.:

(1)
$$𝐯=𝐮^{q^\theta }𝐮.$$

The condition $`gcd(h,q^n1)=1`$ is equivalent to requiring that the map $`𝐮𝐮^h`$ on $`𝕂`$ is $`\mathit{1}\mathit{1}`$; its inverse is the map $`𝐮𝐮^h^{},`$ where $`h^{}`$ is the inverse multiplicative of $`h`$ modulo $`q^n1`$.

In addition, Alice chooses two secret affine transformations, i.e., two invertible matrices $`A=\{A_{ij}\}`$ and $`B=\{B_{ij}\}`$ with entries in $`𝔽_q`$, and two constant vectors $`𝐜=(c_1,\mathrm{},c_n)`$ and $`𝐝=(d_1,\mathrm{},d_n)`$.

Now she sets:

(2)
$$𝐮=A𝐱+𝐜and𝐯=B𝐲+𝐝.$$

Recall that the operation of raising to the $`q^k`$-th power in $`𝕂`$ is an $`𝔽_q`$-linear transformation. Let $`P^{(k)}=\{p_{ij}^{(k)}\}`$ be the matrix of this linear transformation in the basis $`\beta _1,\beta _2,\mathrm{},\beta _n`$, i.e.:

(3)
$$\beta _i^{q^k}=\underset{j=1}{\overset{n}{}}p_{ij}^{(k)}\beta _j,p_{ij}^{(k)}𝔽_q,$$

for $`1i,kn`$. Alice also writes all products of basis elements in terms of the basis, i.e.:

(4)
$$\beta _i\beta _j=\underset{\mathrm{}=1}{\overset{n}{}}m_{ij\mathrm{}}\beta _{\mathrm{}},m_{ij\mathrm{}}𝔽_q,$$

for each $`1i,jn`$. Now she expands the equation (1). So she obtains a system of equations, explicit in the $`v`$, and quadratic in the $`u`$. She uses now her affine relations (2) to replace the $`u,v`$ by the $`x,y`$. So she obtains $`n`$ equations, linear in the $`y`$, and of degree $`2`$ in the $`x`$. Using linear algebra, she can get $`n`$ explicit equations, one for each $`y`$ as polynomials of degree $`2`$ in the $`x`$.

Alice makes these equations public. Bob to send her a message $`(x_1,x_2,\mathrm{},x_n)`$, substitutes it into the public equations. So he obtains a linear system of equations in the $`y`$. He solves it, and sends $`𝐲=(y_1,y_2,\mathrm{},y_n)`$ to Alice.

To eavesdrop, Eve has to substitute $`(y_1,y_2,\mathrm{},y_n)`$ into the public equations, and solve the nonlinear system of equations for the unknowns $`x`$.

When Alice receives y, she decrypts:

$`y_1,y_2,\mathrm{},y_n`$
$``$
$`𝐯=B𝐲+𝐝`$
$``$
$`𝐯={\displaystyle v_i\beta _i}`$
$``$
$`𝐮=𝐯^h^{}`$
$``$
$`𝐱=A^1(𝐮𝐜).`$
In Eurocrypt $`{}_{}{}^{}88`$ \[imai2\], Imai and Matsumoto proposed a digital signature algorithm for their cryptosystem. At Crypto $`{}_{}{}^{}95`$, Jacques Patarin \[Patarin95\] showed how to break this cryptosystem. He noticed that if one takes the equation $`𝐯=𝐮^{q^\theta +1}`$, raises both sides on the $`(q^\theta 1)`$-th power, and multiplies both sides by $`\mathrm{𝐮𝐯}`$, he gets the equation $`\mathrm{𝐮𝐯}^{q^\theta }=𝐮^{q^{2\theta }}𝐯`$ that leads to equations in the $`x`$, $`y`$, linear in both sets of variables. Essentially the equations do not suffice to identify uniquely the message, but now even an exhaustive search will be feasible. The system was definitively insecure and breakable, but its ideas inspired a whole class of public key cryptosystems and digital signatures based on structural identities for finite field operations \[HFE, moh, koblitz, Patarin96, patarin96hidden, gou-pat1\].

Actually, the security of this class lies on the difficulty of the problem of solving systems of polynomial equations. This problem is hard iff the equations are randomly chosen. All manipulations aim to make equations seem like that. If they really were random, the problem is hard to Alice, too.

Our paper is organized as follows. In the next section we develop our own, new cryptosystem. Alice builds her public key by manipulations as above, starting from a certain bivariate polynomial. All of Alice’s manipulations are meant to hide from Eve this polynomial. It is the most important part of the private key. Its knowledge reduces decryption to the practically easy problem of solving a single univariate polynomial.

In the third we discuss some security issues. There we explain that practically all bivariate nonlinear polynomials are good to us to give raise to a public key. This plentitude of choices is an important security parameter.

In the fourth section we provide our cryptosystem with a digital signature algorithm. In the fifth one we provide one more encryption protocol, now a probabilistic one, in the sense that to the same cleartext correspond zero, one, or more cyphertexts.

In the sixth one we discuss some more variations. Essentially, we replace the single bivariate polynomial by an ideal of a small size.

In the seventh section we mention what Shannon \[stinson\] calls Unconditionally Secure Cryptosystems. Actually, this class of cryptosystems is considered an exclusive domain of private key cryptography. This is due mostly to the unhappy state of art of public key cryptography.

In the eighth one we extend our constructions to differential fields of positive characteristic. We hope they are the suitable environment for unconditionally secure public key cryptosystems.

## 2. A New Cryptosystem

### 2.1. Key Generation

Alice chooses two finite fields $`𝔽_q<𝕂`$, and a basis $`\beta _1,\beta _2,\mathrm{},\beta _n`$ of $`𝕂`$ as an $`𝔽_q`$-vector space. Next she takes a generic (for now) randomly chosen bivariate polynomial:

(5)
$$f(X,Y)=\underset{ij}{}𝐚_{ij}X^iY^j$$

in $`𝕂[X,Y]`$, such that she is able to find all its roots in $`𝕂`$ with respect to $`X`$; $``$ $`Y𝕂`$, if any. For the range of $`i`$ employed, this is nowadays considered a relatively easy problem. Further, $`f(X,Y)`$ is subject to other few constraints, that we make clear at the opportune moment.

In transforming cleartext into ciphertext message, Alice will work with two intermediate vectors, $`𝐮=(u_1,\mathrm{},u_n)`$ and $`𝐯=(v_1,\mathrm{},v_n)`$; $`𝐮,𝐯𝕂`$. She sets:

(6)
$$\underset{ij}{}𝐚_{ij}𝐮^i𝐯^j=0.$$

For $`𝐚_{ij}0`$, she sets somehow:

(7)
$$i=\underset{k=1}{\overset{n_i}{}}q^{\theta _{ik}},j=\underset{k=1}{\overset{n_j}{}}q^{\theta _{jk}},$$

where $`\theta _{ik},\theta _{jk}n_i,n_j,_{}`$. Here somehow means that (7) need not be the $`q`$-ary representation of $`i`$, $`j`$. Indeed, there is no reason for it to be. We allow to each $`i`$ both opportunities: to be or not to be. Doing so we increase our choices, whence the random-looking of the public key. In any fashion, what we are dealing with, are nothing but identities.

Next Alice substitutes the (7) to the exponents in (6), obtaining:

(8)
$$\underset{ij}{}(𝐚_{ij}exp(𝐮,\underset{k=1}{\overset{n_i}{}}q^{\theta _{ik}})exp(𝐯,\underset{k=1}{\overset{n_0}{}}q^{\theta _{jk}}))=0;$$

that is:

(9)
$$\underset{ij}{}(𝐚_{ij}\underset{k=1}{\overset{n_i}{}}𝐮^{q^{\theta _{ik}}}\underset{k=1}{\overset{n_j}{}}𝐯^{q^{\theta _{jk}}})=0.$$

Recall that the operation of raising to the $`q^k`$-th power in $`𝕂`$ is an $`𝔽_q`$-linear transformation. Let $`P^{(k)}=\{p_\mathrm{}m^{(k)}\}`$ be the matrix of this linear transformation in the basis $`\beta _1,\beta _2,\mathrm{},\beta _n`$, i.e.:

(10)
$$\beta _i^{q^k}=\underset{j=1}{\overset{n}{}}p_{ij}^{(k)}\beta _j,p_{ij}^{(k)}𝔽_q;$$

for $`1i,jn`$. Alice also writes all products of basis elements in terms of the basis, i.e.:

(11)
$$\beta _i\beta _j=\underset{k=1}{\overset{n}{}}m_{ijk}\beta _k,m_{ijk}𝔽_q;$$

for $`1i,jn`$.

Now she substitutes $`𝐮=(u_1,u_2,\mathrm{},u_n)`$, $`𝐚_{ij}=(a_{ij1},a_{ij2},\mathrm{},a_{ijn})`$, $`𝐯=(v_1,v_2,\mathrm{},v_n)`$, and the identities (10), (11) to (9), and expands. So she obtains a system of $`n`$ equations of degree $`t`$ in the $`u`$, $`v`$, where:

(12)
$$t=max\{n_i+n_j:𝐚_{ij}0\}.$$

Every term under $`\mathrm{\Sigma }`$ in (7) contributes by one to the degree in the $`u`$ of the polynomials.

Here we pause to give some constraints on the range of $`i`$, $`j`$ in (6). The aim of this section is to generate a set of polynomials; linear in a set of variables, and nonlinear in another one. For that purpose, we relate (6) and (7): $`𝐚_{ij}0`$ $`\{n_i>1`$, $`n_j=1\}`$.

On the other side, the size of public key will be $`𝒪((2n)^{t+1})`$. So, it grows polynomially with $`n`$, and exponentially with $`t`$. Therefore, we are interested to keep $`t`$ rather modest, e.g., $`t=2,3`$ or so. So, we have to choose $`i`$, $`j`$ in (5), (7) in order to keep $`t`$ under a forefixed bound.

Next, Alice chooses $`A=\{A_{ij}\},B=\{B_{ij}\}GL(𝔽_q)`$, $`𝐜,𝐝𝕂`$, and sets:

(13)
$$𝐮=A𝐱+𝐜,𝐯=B𝐲+𝐝,$$

where $`𝐱=(x_1,x_2,\mathrm{},x_n)`$, $`𝐲=(y_1,y_2,\mathrm{},y_n)`$ are vectors of variables.

Now she substitutes (13) to the equations in the $`u`$, $`v`$ above, and expands. So she obtains a system of $`n`$ equations of degree $`t`$ in the $`x`$, $`y`$; linear in the $`y`$, and nonlinear in the $`x`$.

After the affine transformation, in each equation appear terms of each degree, from zero to $`t`$; before not. This is its use; to shuffle terms coming from different monomials of (9).

At this point, we are ready to define the cryptosystem.

### 2.2. The Protocol

With the notations adopted above, we define the HPE Cryptosystem (Hidden Polynomial Equations) as the public key cryptosystem such that:

* The public key is:

+ The set of the polynomial equations in the $`x`$, $`y`$ as above;
    + The field $`𝔽_q`$;
    + The alphabet: a set of elements of $`𝔽_q`$.
* The private key is:

+ The polynomial (5);
    + $`A`$, $`B`$, $`𝐜`$, $`𝐝`$ as in (13);
    + The identities (6) to (11);
    + The field $`𝕂`$.
* Encryption:

Bob separates the cleartext $`M`$ by every $`n`$ letters. If needed, he completes the last string with empty spaces. Next he takes an $`n`$-tuple $`𝐱=(x_1,x_2,\mathrm{},x_n)`$ of $`M`$, substitutes it to the $`x`$ in the public equations, solves with respect to the $`y`$, and sends $`𝐲=(y_1,y_2,\mathrm{},y_n)`$ to Alice. We assume here that the solutions exist, and postpone the case they do not.
* Decryption:

Alice substitutes $`𝐯=B^1(𝐲𝐝)𝕂>𝔽_q`$ in (6), and finds all solutions within $`𝕂`$. There is at least one. Indeed, if $`𝐱`$ is Bob’s cleartext, $`𝐮`$ as in (13) is one. For each solution $`𝐮`$, she solves:

(14)
$$𝐱=A^1(𝐮𝐜),$$

and represents all solutions in the basis $`\beta _1,\beta _2,\mathrm{},\beta _n`$. It takes a Chinese Remainder Theorem. With probability $`1`$, all results but one, Bob’s $`(x_1,x_2,\mathrm{},x_n)`$, are gibberish, or even stretch out of the alphabet.

### 2.3. Remarks

#### 2.3.1.

The risc of uncertain decryption is quite virtual. It equals the probability that another sensate combination of letters $`𝐱`$ satisfies (14) for any root $`𝐮`$ of (6) for Bob’s $`𝐲`$, besides the good one that always does. Afterwards, the undesired solution has to join well with the other parts of the decrypted message.

#### 2.3.2.

The main suspended question is that of existence of solutions. Well, Bob succeeds to encrypt a certain message x iff Alice’s equation (6) has solutions for u as in (13) for that x. Alice’s polynomial is a random one. It is a well-known fact from algebra that the probability that a random polynomial of degree $`m`$ with coefficients upon a field $`𝔽_{q^n}`$ has a root in it is about $`1\frac{1}{e}63.2\%`$ \[koblitz, marcus\]. Now the remedy is probabilistic. Alice renders the alphabet public with letters being sets of $`𝔽_q`$. Bob writes down a plaintext and gives start to encryption. If he fails, he substitutes a letter of the cleartext with another one of the same set, and retries.

After $`s`$ trials, the probability he does not succeed is $`\frac{1}{e^s}`$; sufficiently small for the algorithm to be trusted to succeed.

#### 2.3.3.

The other problem is that Alice may have to distinguish the right solution among a great number of them. Here we propose a first remedy. Her number of solution is bounded above by the degree in $`X`$ of $`f`$. So, it is beter to her to keep this degree moderate. Later in this paper in other settings, there will be other remedies, too.

There are no bounds on the degree in $`Y`$. It can be taken whatsoever huge.

#### 2.3.4.

Solving univariate polynomial equations is used by Patarin, too \[patarin96hidden, Wolf:02:Thesis\]. He takes a univariate polynomial:

$$f(x)=\underset{i,j}{}\beta _{ij}x^{q^{\theta _{ij}}+q^{\phi _{ij}}}+\underset{i}{}\alpha _ix^{q^{\xi _i}}+\mu _0,$$

and with manipulations like ours, both the same as Imai-Matsumoto \[imai1\], he gets his public key; a set of quadratic equations. He uses two affine transformations to shuffle the equations. We claim that the first one adds nothing to the security.

The bigger the degree of $`f`$ is, the more the public key resembles a randomly chosen set of quadratic equations. So, it is a security parameter. On the other side, it slows down decryption, principally by adding a lot of undesired solutions. To face that second problem, to the public key are added other, randomly chosen, equations. This is its Achilles’ heel. It makes the public key overdefined, therefore subject to certain facilities to solve \[ckps\]. So, it weakens the trapdoor problem.

We do not add equations to discard undesired solutions. So, we are not subject to overdefined stuff. If in certain variations we do add, we need to add less equations, however. We label wrong solutions those that after decrypted do not make sense, or stretch out of the alphabet.

Afterall, all decrypted texts will howsoever be in a comprehensible language (to someone or some wedget). As $`n`$ grows, it is less possible to have more than one meaningful solution. Besides, any monkey solution that appears to Alice, appears to Eve, too. Furthermore, Eve may have more meaningful solutions. If desired, other tests can be introduced for that purpose. There is no need, however. The solutions, the good one and the bad ones, are very few; no more than $`m`$.

A big advantage of our settings is that we need a lower degree polynomial in $`X`$. So, we make the presence of undesired solutions virtual. Decryption is a pure linear algebra matter.

What is most important, we have now a practically infinite range of choices of $`f`$. This is not Patarin’s case. There the choices are bounded below because of being easy to attack cases, and above because of being impractical to legitimate users.

The only few constraints we put on its monomials aim to:

* keep public key equations linear in the $`y`$;
* have less undesired solutions in decryption process;
* keep the size of public key moderate;
* keep all public key equations nonlinear in th $`x`$.

We can take the degree in y unreasonably high. It gives no trouble to us. It suffices that all the powers of y that appear in the monomials of $`f`$ are powers of $`q`$, so the public equations come linear with respect to the $`y`$.

A new facility now is that we can take lower degree in x, as multiple linear attack does not anymore apply, hopingly.

The constraint that all public key equations must be nonlinear in the $`x`$ is the only non-negotiable one. Indeed, if Alice violates it, the trapdoor problem becomes fatally easy to Gröbner techniques.

Back to the degree in the $`y`$ of the public key. Assume that the public key equations are not linear in the $`y`$. Once Bob substitutes the $`x`$ in the public equations, he now is not challenged to solve a nonlinear system of equations. He is only required to find one solution of it. This can be done within polynomial time with respect to the total degree of the system. Later we give settings to keep public key nonlinear of modest degree in the $`y`$.

Each of such solutions (if any) is encryption to the same cleartext. So we have set up a probabilistic encryption protocol. To a single cleartext may correspond zero, one, or more ciphertexts.

So, in conclusion, Alice is allowed to take for the construction of her public key any damned bivariate polynomial. Indeed, we later argue that $`f`$ can quite well be a multivariate polynomial.

We hope this plentitude of choices is a spoil-sport to Eve.

## 3. Security Issues

Apparently, the only things Eve knows, are the system of public equations, and the order of extension. By brute force, she has to take $`(y_1,y_2,\mathrm{},y_n)`$, to substitute it in the public key equations, to solve in $``$, or maybe $`[\alpha ]`$, and to take the sensate solution. Almost surely, there is only one good solution among those that she finds. She has to find it among $`t^n`$ of them. However, the main difficulty to her is just solving the system. Supposedly, it will pass through the complete computation of Gröbner basis. It is a well-known hard problem. The complexity of computations upon a field grows at most twice exponentially with respect to the number of variables, and in the average case, exponentially.

So, it is better to take $`n`$ huge. This diminishes the probability that Alice confuses decryption, however close to zero, and, what is most important, this renders Eve’s task harder.

Alice and Bob will have to solve sets of bigger systems of linear equations, and face Chinese Remainder Theorem for bigger $`n`$.

There exist well-known facilities \[ckps\] to solve overdefined systems of equations. Unlike most of the rest, our public key is irrendundant, so it is not subject to such facilities.

Now, by exhaustive search we mean that Eve substitutes the $`𝐲`$ in the public equations, and tries to solve it by substituting values to $`𝐱`$.

If we have $`d`$ letters each of them being represented by a single element of $`𝔽_q`$, the complexity of an exhaustive search is $`𝒪(d^n)`$. It is easy for Alice to render exhaustive search more cumbersome than Gröbner attack. The last one seems to be the only choice to Eve.

We did not find any Known Cleartext Attack to our cryptosystem.

Eve may engineer cleartext$``$ciphertext analyses, seeking for invariants or regularities there, helpful for an attack \[patarin96hidden\]. All the identities we use, mean totousle any such regularity, and to disguise from Eve any hint on $`i`$, $`j`$, and on the entries of $`A`$, $`B`$, $`𝐜`$, $`𝐝`$, and the $`𝐚_{ij}`$; that she may use for such an attack.

The complexity of the trapdoor problem is $`𝒪(t^n)`$, the size of public key $`𝒪(n^{t+1})`$. This fully suggests the values of parameters. $`n=100`$, $`t=2,3,4`$ would be quite good choices.

Obviously, infinitely many bivariate polynomials give raise to the same public key. Indeed, fixed the ground field, the degree of extension $`n`$, and the degree of public key equations, we have a finite number of public keys. On the other hand, there are infinitely many bivariate polynomials that can be used like private keys.

On how does it happen, nothing is known. If ever found, any such regularity will only weaken the trapdoor problem.

## 4. A Digital Signature Algorithm

Assume that we are publicly given a set of hash functions that send cleartexts to strings of integers of fixed length $`n_B`$. For the only purpose of signing messages for Alice, Bob builds a cryptosystem as above with $`q_B`$ prime, and $`[𝕂_B:𝔽_{q_B}]=n_B`$. He to sign a message $`M`$:

* calculates $`H(M)=(y_1,y_2,\mathrm{},y_{n_B})=𝐲𝕂_B`$;
* finds one solution (if any; otherwise, see section 2.3.2.) u of $`f_B(𝐮)=𝐲`$ in $`𝕂_B`$.
* calculates $`𝐱=B^1(𝐮𝐜_B)`$;
* appends $`𝐱=(x_1,x_2,\mathrm{},x_{n_B})`$ to $`M`$, encrypts, and sends it to Alice. $`(x_1,x_2,\mathrm{},x_{n_B})`$ is a signature to $`M`$.

To authenticate, Alice first decrypts, then she:

* calculates $`H(M)=(y_1,y_2,\mathrm{},y_{n_B})`$;
* substitutes $`(x_1,x_2,\mathrm{},x_{n_B})`$, $`(y_1,y_2,\mathrm{},y_{n_B})`$ to Bob’s public equations;
* so she gets an $`n_B`$-tuple of integers. If they all reduce to zero modulo $`q_B`$, she accepts the message; otherwise she knows that Eve has been causing trouble.

If Eve tries to impersonate Bob and send to Alice her own message with hash value $`𝐲=(y_1,y_2,\mathrm{},y_{n_B})`$, then to find a signature $`(x_1,x_2,\mathrm{},x_{n_B})`$, she may try to find one solution of Bob’s system of equations for y. We trust on the hardness of this problem for the security of authentication.

## 5. A Probabilistic Encryption Protocol

With the ideas described above, we are going to set up now a probabilistic protocol such that only the legitimate users can send messages to which-another. Mean, the message is meaningful iff there are no intruders. Its being meaningful is the signature itself.

Here is the shortest possible description. Let $`F_A`$ and $`F_B`$ be Alice’s and Bob’s public keys functions respectively, where $`n_A=n_B`$. To send a message x to Alice, Bob sends her a random (this randomness is the probabilistic pattern) element of $`F_A(F_B^1(𝐱))`$, that she can decrypt by calculating $`F_B(F_A^1(F_A(F_B^1(𝐱))))`$. So if $`F_A(F_B^1(𝐱))\mathrm{}`$. Otherwise, the approach is probabilistic, as in the previous section.

Here is the extended description. Each (English, e.g.) letter (or some of them, only) is represented by a set of few (two, e.g.) elements of the field, or strings of them. For ease of explanation, Bob’s public equations are linear in the $`x`$, and of higher degree in the $`z`$.

Bob writes down the cleartext $`𝐱`$ and finds one solution of:

(15)
$$𝐱=𝐛_r𝐳^r+𝐛_{r1}𝐳^{r1}+\mathrm{}+𝐛_0.$$

If there are no solutions, Bob changes a representant of a letter, and retries. Probability issues are discussed in the previous section.

Now Bob takes the solution z of (15), and applies:

(16)
$$𝐲^{}=B^1(𝐳𝐜_B).$$

Next he takes $`𝐲^{}`$, substitutes in Alice’s public equations. So he obtains a tuple y, that he sends to Alice. This is the ciphertext.

Each of other solutions of (15) give raise to other encryptions of the same cleartext.

Alice now to decrypt, solves her equation for y within her field $`𝕂`$. There is at least one solution. Next she applies her inverse affine transformation to all (few) solutions, and substitutes them all on Bob’s public equations. Of that procedure all, Alice now discards all meaningless solutions, and takes the meaningful one.

What is the trapdoor problem now? Well, on authentication matter, nothing new. Eve has the same chances to forge here that she had before. Recall that this kin of signatures is already best with respect to the other ones.

On security, instead, there is a very good improvement. By brute force, Eve has to take the ciphertext, substitute on Alice’s public key, find all solutions, and substitute them all on Bob’s public key; then take the sensate ones. This is worse than exhaustive search of previous cryptosystems.

Now, what does here really mean exhaustive search? Eve now has to search through all the elements of the common public ground field, not just through all the alphabet. So, opting for this protocol, we can put a lot of constraints on alphabet, in order to discard far easier the undesired solutions, without rendering the public key overdefined.

She sets up such $`n`$-tuples, checks whether they are solutions of Alice’s public key for Bob’s ciphertext y substituted to the variables $`y`$. If yes, she substitutes to Bob’s public key, and checks whether does it make sense.

What can linear multiple attack or quadratic attack \[patarin96hidden\] do in these new settings?

Apart all, we save space and calculi. We do not need any more the calculi and space of signature.

This protocol can be used for multiple encryption, too.

Let us suppose that the letters are strings of a fixed length. Well, here Alice can impose that not all strings are letters. So, in decryption she discards a priori the solutions that contain non-letters. Doing so, she actually has a single good solution of her polynomial, and saves herself the effort of appealing to other tricks. In all the other schemes throughout, such a trick fatally weakens the exhaustive search.

## 6. Hidden Ideal Equations

Instead of a single bivariate polynomial, Alice may choose to employ an ideal of a very modest size. She separates the variables she employs into two sets, $`\{X_i\}`$, $`\{Y_j\}`$; one for encryption, one for decryption. She may decide to leave one of the equations employed of higher degree in the $`\{Y_j\}`$ after manipulations, so she gives raise to a probabilistic encryption protocol. Alice’s parameters are:

* $`n=[𝕂:𝔽_q]`$;
* the number $`s_1`$, $`s_2`$ of variables $`\{X_i\}`$, $`\{Y_j\}`$, respectively;
* the number $`r`$ of private equations.

So, the number of public key equations is $`nr`$. The number of the variables $`x_{ij}`$ is $`ns_1`$, and that of the $`y_{kl}`$ is $`ns_2`$.

Alice’s number of variables, the $`\{X_i\}`$, is insignificant so far, so she is supposed to be able to appeal to Gröbner stuff in order to solve her system of equations within the field of coefficients for Bob’s $`\{Y_j\}`$.

What is most important here and throughout, if Bob succeeds to encrypt, Alice does always succeed to decrypt.

For ease of treatment, assume now that Alice does not apply affine transformations to her variables. Bob fails encryption for a certain cleartext $`(X_1,\mathrm{}X_{s_1})`$ iff Alice’s private ideal has no solutions in the $`Y`$ for such an $`(X_1,\mathrm{}X_{s_1})`$. Alice’s private ideal is a random one. If she takes $`rs_2`$, the probability that it has no solutions is $`0`$, and $`1`$ for $`r>s_2`$. So, it suffices that Alice takes $`rs_2`$. The critical cases that may supervene are faced simply changing alphabet.

With slight changes, this reasoning holds in the case that Alice applies affine transformations, too.

The real problem is indeed that the solutions to Alice may be too many; and in any case finitely many, as the base field is finite. The best remedy to that is that Alice takes $`r=s_1`$. So, the ideal that she obtains after substitution of Bob’s ciphertext is zerodimensional (quite easy to cause it happen), and the number of solutions is bounded above by the total degree of the system. So, she can contain the number of solutions by taking the total degree in the $`\{X_i\}`$ modest, and however each of them nonlinear.

Alice can take all equations of very low degree in the $`X`$, and then transform that basis of the ideal they generate to another one of very high degrees in the $`X`$. So she has a low Bezout number of the ideal, and higher degrees in the $`X`$, and transformations as above can take place. If she takes the first basis linear, the number of solutions of her equations reduce to one: Bob’s cleartext.

As soon as $`r>s_1`$, the public key becomes overdefined.

Alice applies a permutation to the equations and a renumeration to the variables before publishing her key, so Eve does not know how are they related. She may apply affine transformations, or may not, or may apply to only some of the $`X_i`$, $`Y_j`$; at her discretion.

If $`s_1<s_2`$, the size of the ciphertext is bigger than that of cleartext, and nothing else wrong. By this case, encryption is practically always probabilistic. Indeed, even when the equations are linear with respect to the $`y_{kl}`$, since there are more variables than equations, the solutions exist, and are not unique.

Actually, Alice can take $`s_2`$ rather huge. She may choose to manipulate some of the $`Y_j`$ within a subfield of $`𝕂`$, rather than within $`𝕂`$. Doing so, she allows herself a big $`s_2`$, and a contained size of the ciphertext. The number of the variables $`y_{kl}`$ now is no more $`ns_2`$.

### 6.1.

Now the size of the public key is $`𝒪(s_1(n)^{t+1})`$, and the complexity of the trapdoor problem is $`𝒪(t^{ns_1})`$.

It is true that throughout the size of public key grows polynomially with $`n`$, but before $`n`$ becomes interesting, the public key is already quite cumbersome. So, opting for the choices of this section we have reasonable security with much smaller values of $`n`$. $`n=20`$, or so, actually are quite good. We are allowed some more values of $`t`$, too.

### 6.2.

There exist classes of ideals called with doubly exponential ideal membership property \[swanson\]. These are the ideals for which the calculus of a Gröbner basis cannot be done within exponential time on the number of variables, i.e., it can be done within doubly exponential time on the number of variables. It is very interesting to know whether can we employ them in some fashion in this class of cryptosystems. In any fashion, this is the theoretical limit for employing solving of polynomial systems of equations in public key cryptography.

## 7. Some Considerations

The idea of public key cryptography was first proposed by Diffie and Hellman \[pkc\]. Since then, it has seen several vicissitudes \[odlyzko\].

A trapdoor function is a map from cleartext units to ciphertext units that can be feasibly computed by anyone having the public key, but whose inverse function cannot be computed without knowledge of the private key:

* either because (at present, publicly) there is no theory to do it;
* or the theory exists, but the amount of calculations is deterring.

Cryptosystems with trapdoor problems of the first kin are what Shannon \[stinson\] calls Unconditionally Secure Cryptosystems.

Actually, the aim is to make trapdoor problems be equivalent to time-honoured hard mathematical problems. However, being of a problem hard or undecidable implies nothing about the security of the cryptosystem \[barkee, odlyzko\]. Recall that of all schemes ever invented, only two of them, $`RSA`$ \[rsa\] and ECDL \[koblitz\], are going to be broken (or, at least, are going to become impractical) by solving the hard problems they lie upon. The rest of them have been broken with theories of no use to solve their hard problem. So, once more, it may happen to be proved that solving systems of differential&integral equations is undecidable, nevertheless several cryptosystems built upon them may be easy to break rather than secure.

The author is very fond of the idea of public key cryptography, and believes howsoever in new developments that will make it fully suffice for all purposes.

Actually, one tendency is that of investigating poor structures, mean, structures with less operations, like groups, semigroups with cryptosystems upon the word problem \[anshel, yamamura, hughes\]. Yamamura’s paper \[yamamura\] can be considered pioneering on secure schemes. Unfortunately, its scheme is still uneffective.

William Sit and the author are investigating cryptosystems upon other algebraic structures. We are investigating among other things whether is it possible to build effective secure schemes upon differential fields of positive characteristic. We hope that cryptography will arouse new interests on differential and universal algebra, too, as it did in number theory and arithmetic geometry. One reason of optimism is that in universal algebra one can go on further with new structures and hard or undecidable problems forever. Until now we have appealed to only the unary and binary arithmetic operations.

## 8. Generalizations on Differential Fields

Differential algebra is born principally due to the efforts of Ritt \[ritt\] to handle differential equations by means of algebra. Actually, a differential field is a field with a set of unary operations  called derivatives that replace an element of the field with another one such that $`(a+b)^{}=a^{}+b^{}`$ and $`(ab)^{}=ab^{}+a^{}b`$.

Good references in the topic are \[kolchin, sit2, ritt, sadik, kaplansky\]. Kaplansky’s book is probably the best introduction in the topic.

It is possible222Most of considerations given in this section are suggestions of professor Sit through private communications. to generalize the schemes given throughout using differential polynomials instead of (5). Take $`𝕂`$ to be a finite differential field extension of a differential field $`𝔽`$ of positive characteristic333In zero characteristic numerical analysis tools seriously affect security, or at least constrain us to more careful choices. We shall not dwell on this topic here.. Any such $`𝕂`$ is defined by a system of linear homogeneous differential equations, and there are structural constants defining the operations for the derivations (one matrix for each derivation), as well for multiplication.

One can now replace (5) with a differential polynomial. The scheme works verbatim. One can take (5) to be of higher order and degree, that is ok too, just like the algebraic case. Euler, Clairaut, or any of other well-studied classes of equations, or their compositions; each of them fully suffice.

The techniques described above for polynomials, if applied to differential polynomials, will definitely make it much harder to attack any protocol developed. Any affine transformation (by this is meant a linear combination of the differential indeterminates with not-necessarily constant coefficients, and this linear combination is then substituted differentially in place of the differential indeterminates) will not only even out the degrees, but also the orders of the various partials, and making the resulting differential polynomial very dense.

However, there is one thing to caution about: any time one specifies these structural matrices, they have to satisfy compatibility equations. In the algebraic case, it is the relations between $`P^k=\{p_{ij}^{}{}_{}{}^{(k)}\}`$ in (10) and $`M_{\mathrm{}}=\{m_{ij\mathrm{}}\}`$ in (11). The $`P^k`$ are simply determined uniquely by $`M_{\mathrm{}}`$, given the choices implicitely defined in (11).

It is very interesting to know in the algebraic case whether the system of equations Alice obtains is invariant under a change of basis, all other settings being equal. There is probably some group of matrices in $`GL(n,q)`$ that can do that. Such a knowledge may be used to build attacks to all schemes of $`HFE`$ class.

In the differential case there is a similar action called Loewy action, or the gauge transformation. For ordinary differential equations, two matrices $`A`$, $`B`$ are Loewy similar if there is an invertible matrix $`K`$ such that $`A=\delta KK^1+KBK^1`$. Using this action, one can classify the different differential vector space structures of a finite dimensional vector space. There is also a cyclic vector algorithm to find a special basis, so that the differential linear system defining the vector space becomes equivalent to a single linear $`ODE`$.

If no other problems arise for the differential algebraic schemes, there is however one caution more for them to be unconditionally secure. We have to avoid the exhaustive search. For that, Alice has to publish a finite alphabet where each letter is represented by an infinite set, disjoint sets for different letters. This is possible in differential fields, as they are infinite. Alice renders the sets public parametrically, as differential algebraic functions of elements of the base differential field, and parameters, e.g., in $``$. Bob chooses a letter, gives random values to parameters, obtains one representant of the letter, and proceeds as above. In any case, if $`\mu `$ is the order of public equations, any two elements $`\mathrm{\Xi }`$, $`\mathrm{\Theta }𝔽`$ such that $`(\mathrm{\Xi }\mathrm{\Theta })^{(\mu )}=0`$ must represent the same letter, if any.

The main care for Alice is that the public key equations must not fall into tractable classes by well-known means, such as linear algebra.

In the algebraic case such constructions do not make sense. Eve can anyway appeal to Gröbner attack. Besides, in any fashion such data enable her to guess $`q`$.

The size of the public key now is actually $`𝒪(n^{to+1})`$, where $`o`$ is the order of public key equations. Quite explosive. However, a first tool to contain it is the low characteristic of the field. So, we see a lot of monomials reduce to zero. The best consolation is that we do not have to go far away with parameters. The trapdoor problem is simply undecidable. $`n=20`$ would fully suffice. Such a value is needed more in order to avoid uncertain decryption, however less probable in differential fields, as the range of solutions is infinite, than for growing security. Besides, if there was found some attack for the $`HDPE`$ (Hidden Differential Polynomial Equations) scheme, it will work better with $`HPE`$. As of now, $`HDPE`$ trapdoor problem seems undecidable, and the scheme effective. The author is working to come up with concrete examples of this kind of cryptosystems. Unfortunately, everything in the topic is still handmade, and therefore rather time-consuming.

### Acknowledgments.

The author would like to thank Don Coppersmith, Patrizia Gianni, Teo Mora, Massimiliano Sala, and Barry Trager for many suggestions and fruitful discussions. The author is particularly indebted to William Sit for several comments and improvements on earlier drafts, and to his advisor, Carlo Traverso.