text:
# Zeros of systems of exponential sums and trigonometric polynomials

## 1. Motivation and Summary

Algebraic geometry is concerned with the study of zero sets of algebraic polynomials $`f(x)=_\alpha c_\alpha x^\alpha `$, where $`\alpha =(\alpha _1,\mathrm{},\alpha _n)^n`$, $`x=(x_1,\mathrm{},x_n)`$. If we allow the exponents $`\alpha `$ to be real vectors, $`f`$ becomes a multi-valued function. This can be remedied by an exponential change of variables $`x_i=\mathrm{exp}2\pi z_i`$, after which we obtain a single-valued exponential sum of the form $`_\alpha c_\alpha \mathrm{exp}2\pi \alpha z`$, where $`\alpha =(\alpha _1,\mathrm{},\alpha _n)^n`$ is the vector of frequencies, $`z=(z_1,\mathrm{},z_n)`$ is the vector of variables, and $`\alpha z`$ is the standard scalar product.

Some results from algebraic geometry can be generalized to this wider class of functions. For example, Bernstein’s theorem states that the number of zeros of a generic system of $`n`$ algebraic equations in $`(^\times )^n`$ with a fixed collection of Newton polyhedra is equal to $`n!`$ times the mixed volume of the Newton polyhedra of the system. A system of $`n`$ exponential equations in $`n`$ variables usually has infinitely many isolated zeros, thus one has to study the distribution of these zeros in order to obtain finite invariants. O. Gelfond proved that the mean number of complex zeros of a system of $`n`$ exponential sums in $`n`$ variables with real frequencies whose Newton polyhedra have sufficiently general mutual positions (so-called developed collection of Newton polyhedra), is equal to $`n!`$ times the mixed volume of the Newton polyhedra of the system .

O. Gelfond and A. Khovanskii found a formula for the sum of the values of a Laurent polynomial at the zeros of a system of $`n`$ Laurent polynomials in $`^n`$ with a developed collection of Newton polyhedra . This formula splits into two components. One of them is geometrical and reflects the mutual positions of the Newton polyhedra of the system, while the other component is expressed explicitly in terms of the coefficients of the polynomial and the system.

An exponential change of variables gives a similar formula for exponential sums with rational frequencies. We conjecture that this formula also holds for exponential sums with real frequencies. Here is some evidence for this. If the exponential sum that we are summing up is identically equal to one, the formula follows from combining two results: Gelfond’s generalization of Bernstein’s theorem  and the new formula for the mixed volume . In  the conjectured formula is proved in dimension one.

The conjectured formula, first of all, implies that the mean value exists. We prove the existence not only in the complex situation but also in a very general real setting by providing an integral formula for the mean value.

If the frequencies of the exponential sum that we are summing up are not commensurate with the frequencies of the system the conjectured formula states that the mean value is equal to zero. We show that this is actually true and therefore prove the conjectured formula in most cases. For example, if the exponential sum that we are summing up is a single exponent $`\mathrm{exp}2\pi \alpha z`$, $`z^n,\alpha ^n`$, then the formula is proved for all values of $`\alpha `$ except for a countable set in $`^n`$.

Our arguments represent a combination of real analytic geometry and ergodic theory. They are based on two theorems from completely different parts of mathematics: the cell-decomposition theorem for subanalytic sets (Appendix B) and Weyl’s equidistribution law for multidimensional trajectories in the real torus (Appendix A). Weyl’s equidistribution law for one-dimensional trajectories is a classical theorem which was published in 1916 (see ). The proof of the corresponding law for multidimensional trajectories in  is a direct generalization of Weyl’s original argument.

Remarks. Similar ideas first appeared in Weyl’s papers  and  where he solves the mean motion problem. S. Gusein-Zade and A. Esterov were dealing with related questions and were using close techniques in  to prove the existence of the mean Euler characteristic and mean Betti numbers of level sets and sets of smaller values of a quasiperiodic function. Their motivation is completely different from ours: it comes from the analysis of some models of chaotic behavior appearing in quasicrystal structures. See also a recent paper  of V. I. Arnold for a close discussion.

Acknowledgements. This paper is a part of the author’s Ph.D. thesis  defended in spring 2002. I would like to thank my thesis advisor Askold Khovanskii for stating the problem and for his constant attention to this work.

## 2. The Gelfond–Khovanskii formula and its conjectured generalization to the case of exponential sums with real frequencies

Let $`\mathrm{\Delta }_1,\mathrm{},\mathrm{\Delta }_n`$ be convex polyhedra in $`^n`$ and $`\mathrm{\Delta }=\mathrm{\Delta }_1+\mathrm{}+\mathrm{\Delta }_n`$ be their Minkowski sum.

###### Definition 2.1.

A collection of faces $`\{\mathrm{\Gamma }_j:\mathrm{\Gamma }_j\mathrm{\Delta }_j,j=1\mathrm{}n\}`$ is coordinated if there exists a nonzero linear function on $`^n`$ whose maximum on the polyhedron $`\mathrm{\Delta }_j`$ is attained exactly on the face $`\mathrm{\Gamma }_j\mathrm{\Delta }_j`$ for all $`j=1\mathrm{}n`$. A collection of polyhedra $`\mathrm{\Delta }_1,\mathrm{},\mathrm{\Delta }_n`$ is developed if none of the polyhedra is a vertex, and in any coordinated collection of faces there is at least one vertex.

Fix a system

$$P_1(z)=\mathrm{}=P_n(z)=0,z^n$$

of $`n`$ Laurent equations with a developed collection of Newton polyhedra $`\mathrm{\Delta }_1,\mathrm{},\mathrm{\Delta }_n`$. The Gelfond–Khovanskii formula states that the sum of the values of a Laurent polynomial $`Q`$ over the zeros of the system in $`(^\times )^n`$ is equal to

$$(1)^n\underset{\alpha }{}k_\alpha C_\alpha ,$$

where the summation is performed over the vertices $`\alpha `$ of $`\mathrm{\Delta }=\mathrm{\Delta }_1+\mathrm{}+\mathrm{\Delta }_n`$, $`k_\alpha `$ is the combinatorial coefficient that corresponds to the vertex $`\alpha `$ (the combinatorial coefficients are integers that reflect the mutual position of the Newton polyhedra $`\mathrm{\Delta }_1,\mathrm{},\mathrm{\Delta }_n`$ of the system, see  for the definition), and $`C_\alpha `$ is an explicit Laurent polynomial in the coefficients of $`P=P_1\mathrm{}P_n`$ and $`Q`$. This result was announced in , and a proof was given in .

If we allow the exponents of the Laurent polynomials to be real vectors, we obtain multi-valued functions. After an exponential change of variables they become single-valued exponential sums with real frequencies. Now we make some preparations before we formulate the conjectured generalization of the Gelfond–Khovanskii formula to the case of exponential sums.

Let $`\mathrm{\Lambda }`$ be a finite set in $`^n`$. An exponential sum with the spectrum $`\mathrm{\Lambda }`$ is a function $`F:^n`$ of the form

$$F(z)=\underset{\alpha \mathrm{\Lambda }}{}c_\alpha \mathrm{exp}2\pi \alpha z,$$

where the summation is performed over the frequencies $`\alpha \mathrm{\Lambda }`$, $`z=(z_1,\mathrm{},z_n)^n`$, $`c_\alpha `$ are nonzero complex numbers, and $`\alpha z`$ is the standard scalar product. The Newton polyhedron of an exponential sum is the convex hull $`\mathrm{\Delta }(\mathrm{\Lambda })`$ of its spectrum $`\mathrm{\Lambda }`$.

Fix a system of exponential sums

(2.1)
$$F_1(z)=\mathrm{}=F_n(z)=0,z^n$$

with a developed collection of Newton polyhedra $`\mathrm{\Delta }_1,\mathrm{},\mathrm{\Delta }_n`$ in $`^n`$. Due to the assumption that the collection of the polyhedra is developed, there exists $`R>0`$ such that all the zeros of the system lie in a strip $`S_R\times \text{Im }^n`$, where $`S_R\text{Re }^n`$ is a ball of radius $`R`$ centered at the origin. This implies that all the zeros are isolated (see Theorem 5.1).

Let $`G`$ be another exponential sum with real frequencies. The sum of the values at the zeros of a system in the exponential case is replaced with the result of averaging $`G`$ over the zeros of the system along the imaginary subspace. Let $`\mathrm{\Omega }\text{Im }^n`$ be a bounded set with nonzero volume. For $`\lambda `$, define $`S_\mathrm{\Omega }(\lambda )`$ to be the sum of the values of $`G`$ at the zeros of (2.1) (counting multiplicities) that belong to the strip $`^n\times \lambda \mathrm{\Omega }^n`$.

###### Definition 2.2.

The mean value $`M_\mathrm{\Omega }`$ of $`G`$ over the zeros of the system (2.1) is the limit of $`S_\mathrm{\Omega }(\lambda )/\mathrm{Vol}(\lambda \mathrm{\Omega })`$, as $`\lambda `$ approaches infinity.

The Minkowski sum $`\mathrm{\Delta }=\mathrm{\Delta }_1+\mathrm{}+\mathrm{\Delta }_n`$ is the Newton polyhedron of the product $`F=F_1\mathrm{}F_n`$. The exponent $`\mathrm{exp}2\pi \alpha z`$ that corresponds to a vertex $`\alpha `$ in $`\mathrm{\Delta }`$ appears in $`F`$ with a nonzero coefficient $`d_\alpha `$. Let $`\stackrel{~}{F}=F/(d_\alpha \mathrm{exp}2\pi \alpha z)`$. The constant term of the exponential sum $`\stackrel{~}{F}`$ is equal to one. We define the exponential series for $`1/\stackrel{~}{F}`$ by the formula

$$1/\stackrel{~}{F}=1+(1\stackrel{~}{F})+(1\stackrel{~}{F})^2+\mathrm{}.$$

Since each exponent appears with a nonzero coefficient in a finite number of terms, the coefficients of this series are well-defined. Let $`C_\alpha `$ be the constant term in the formal product of this series and $`(1/d_\alpha )\mathrm{exp}(2\pi \alpha z)Gdet(\frac{F}{z})`$.

###### Conjecture 2.3.

The mean value $`M_\mathrm{\Omega }`$ of an exponential sum $`G`$ over the zeros of the system (2.1) is equal to

$$\frac{1}{(2\pi )^n}\underset{\alpha }{}k_\alpha C_\alpha ,$$

where the summation is performed over the vertices $`\alpha `$ of the Minkowski sum $`\mathrm{\Delta }=\mathrm{\Delta }_1+\mathrm{}+\mathrm{\Delta }_n`$, $`C_\alpha `$ are the constant terms of the series as defined above, and $`k_\alpha `$ are the combinatorial coefficients.

As we mentioned in the summary, this formula is proved in the following three cases:

* The frequencies $`\alpha `$ are rational (the formula is obtained by an exponential change of variables from the algebraic case).
* $`G=1`$ (the formula follows from combining two results: Gelfond’s generalization of Bernstein’s theorem  and the new formula for mixed volume ).
* $`n=1`$ (see ).

In the exponential case even the existence of the mean value is not obvious. We give an integral formula for the mean value which proves the existence not only in the complex situation (Theorem 5.2) but also in a very general real setting (Theorem 4.3).

Assume that the frequencies of the exponential sum that we are summing up are not commensurate with the frequencies of the system. Then there is no constant term in the series defined above, that is, all $`C_\alpha `$ are equal to zero, and the conjectured formula states that the mean value is equal to zero. Using the integral representation, we show in Corollary 5.3 that this is actually true. This proves the conjectured formula in most cases. For example, if the exponential sum that we are summing up is a single exponent $`\mathrm{exp}2\pi \alpha z`$, then we proved the formula for all values of $`\alpha `$ except for a countable set in $`^n`$.

## 3. Averaging over the isolated intersections of a subanalytic set with a dense orbit on $`𝕋^N`$

Here we present a construction and formulate a theorem on which all our results are based. Let $`\stackrel{~}{V}`$ be a subanalytic subset of the real torus $`𝕋^N`$ and $`𝒪`$ a dense multidimensional orbit on $`𝕋^N`$. We prove an integral formula for the mean value of a bounded subanalytic function on $`\stackrel{~}{V}`$ over the isolated intersections of the set $`\stackrel{~}{V}`$ and the orbit $`𝒪`$. This formula, in particular, implies that the mean value always exists.

Let $`x=(x_1,\mathrm{},x_n)^n`$ and $`\phi =(\phi _1,\mathrm{},\phi _N)𝕋^N=^N/^N`$. A linear map $`\mathrm{\Phi }:^n^N`$ defines an action of $`^n`$ on $`𝕋^N`$ by

$$x:\phi (\phi +\mathrm{\Phi }(x))\mathrm{mod}^N.$$

Let $`𝒪`$ be an orbit of this action. We require that the orbits are dense in the torus, which means that there are no integral vectors orthogonal to the $`n`$-plane $`\mathrm{\Phi }(^n)^N`$ (Appendix A). We will also assume that $`\mathrm{\Phi }`$ is injective, that is, the orbits are $`n`$-dimensional.

Let $`\stackrel{~}{V}`$ be a subanalytic subset of $`𝕋^N`$ and $`\stackrel{~}{T}`$ a bounded subanalytic function on $`\stackrel{~}{V}`$. Set $`V:=\mathrm{\Phi }^1(\stackrel{~}{V}𝒪)^𝓃`$, and define a function $`T(x)=\stackrel{~}{T}(\mathrm{\Phi }(x))`$ on $`V`$. Let $`\mathrm{\Omega }`$ be a bounded subset of $`^n`$ with nonzero volume. For $`\lambda >0`$, define $`S_\mathrm{\Omega }(\lambda )`$ to be the sum of the values of $`T`$ at the isolated points of $`V`$ that belong to $`\lambda \mathrm{\Omega }`$.

###### Definition 3.1.

The mean value $`M_\mathrm{\Omega }`$ of $`\stackrel{~}{T}`$ over the isolated intersections of $`\stackrel{~}{V}`$ and $`𝒪`$ (in the topology of $`𝒪`$) is the limit:

$$M_\mathrm{\Omega }=\underset{\lambda \mathrm{}}{lim}\frac{S_\mathrm{\Omega }(\lambda )}{\mathrm{Vol}(\lambda \mathrm{\Omega })}.$$

Let $`p:^NL`$ be the orthogonal projection to the linear subspace $`L^N`$ orthogonal to $`\mathrm{\Phi }(^n)`$. Define $`M_{Nn}(\stackrel{~}{V})`$ to be the smooth $`(Nn)`$-dimensional part of $`\stackrel{~}{V}`$, transversal to the orbit. Here is the precise three step definition (see Figure 1 below where we assume that the direction of the orbit is vertical):

1. Let $`\stackrel{~}{V}_0`$ be the set of all points $`\phi \stackrel{~}{V}`$ such that the intersection of $`\stackrel{~}{V}`$ with the $`n`$-plane through $`\phi `$ parallel to the orbit is locally just the point $`\phi `$ itself.
2. Let $`\stackrel{~}{V}_1`$ be the set of all points $`\phi \stackrel{~}{V}_0`$ such that $`\stackrel{~}{V}_0`$ is a $`C^1`$-manifold of dimension $`Nn`$ in a small neighborhood of $`\phi `$.
3. Let $`M_{Nn}(\stackrel{~}{V})`$ be the set of all points $`\phi \stackrel{~}{V}_1`$ such that the projection $`p`$ is regular at each of these points.

Let $`A^j=\mathrm{\Phi }(e_j)`$, where $`e_j`$ is the $`j`$-th vector in the standard basis for $`^n`$. The quotient of the standard volume form on $`^N`$ by the image under $`\mathrm{\Phi }`$ of the standard volume form on $`^n`$ defines a volume form $`\omega =d\phi (\mathrm{},A^1,\mathrm{},A^n)`$ and an orientation on $`L`$. Then $`M_{Nn}(\stackrel{~}{V})`$ is a $`C^1`$-manifold, and the projection $`p`$ defines an ($`Nn`$)-form $`p^{}\omega `$ on $`M_{Nn}(\stackrel{~}{V})`$.

###### Definition 3.2.

The $`𝒪`$-transversal $`\stackrel{~}{T}`$-weighted volume of $`\stackrel{~}{V}`$ is the integral of the form $`\stackrel{~}{T}p^{}\omega `$ over $`M_{Nn}(\stackrel{~}{V}).`$

###### Definition 3.3.

We say that the fractal dimension of a set $`X`$ is less than or equal to $`n`$ if there exist constants $`c`$ and $`ϵ_0>0`$ such that for any $`0<ϵ<ϵ_0`$, the set $`X`$ can be covered by no more than $`cϵ^n`$ balls of radius $`ϵ`$. We say that the fractal dimension is less than $`n`$ if it is less than or equal to $`n\delta `$ for some positive $`\delta `$.

Our main result is the following theorem:

###### Theorem 3.4.

Let $`𝒪`$ be an $`n`$-dimensional dense orbit in the torus $`𝕋^N=^N/^N`$, $`n<N`$, $`\stackrel{~}{V}`$ a subanalytic subset of $`𝕋^N`$, and $`\stackrel{~}{T}`$ a bounded subanalytic function on $`\stackrel{~}{V}`$. Let $`\mathrm{\Omega }`$ be a bounded subset of $`^n`$ with nonzero volume whose boundary has fractal dimension less than $`n`$.

Then the mean value $`M_\mathrm{\Omega }`$ of $`\stackrel{~}{T}`$ over the isolated intersections of $`\stackrel{~}{V}`$ and $`𝒪`$ always exists and equals the $`𝒪`$-transversal $`\stackrel{~}{T}`$-weighted volume of $`\stackrel{~}{V}`$.

We prove this theorem in Section 6.

## 4. Zeros of systems of trigonometric polynomials

We apply Theorem 3.4 to prove that the mean value of a quasiperiodic trigonometric polynomial $`T`$ over the isolated points of a quasiperiodic semitrigonometric set $`V`$ (a set described by quasiperiodic trigonometric polynomials, see the precise definition below) always exists. We obtain an integral representation for the mean value from which we conclude that in the case when the frequencies of $`T`$ are not commensurate with the frequencies of the quasiperiodic trigonometric polynomials that describe $`V`$, the mean value is equal to zero.

### 4.1. Definitions

Let $`x^n`$. A quasiperiodic trigonometric polynomial is a trigonometric polynomial of the form

$$T(x)=\underset{k=1}{\overset{p}{}}c_k\mathrm{cos}2\pi \alpha _kx+d_k\mathrm{sin}2\pi \alpha _kx$$

with real coefficients $`c_k,d_k`$ and real frequencies $`\alpha _k^n`$. Here $`\alpha _kx`$ is the standard scalar product.

A subset of $`^n`$ is called quasiperiodic semitrigonometric if it is described by quasiperiodic trigonometric equations and inequalities. More precisely, we say that $`V^n`$ is quasiperiodic semitrigonometric if it can be represented in the form:

$$V=\underset{i=1}{\overset{s}{}}\underset{j=1}{\overset{r_i}{}}V_{ij},$$

where each $`V_{ij}`$ is either $`\{x^n|T_{ij}(x)=0\}`$ or $`\{x^n|T_{ij}(x)>0\}`$, for some quasiperiodic trigonometric polynomials $`T_{ij}`$.

Let $`T`$ be a quasiperiodic trigonometric polynomial, $`V`$ a quasiperiodic semitrigonometric set in $`^n`$ described by $`T_{ij}`$, and $`\mathrm{\Omega }`$ a bounded subset of $`^n`$ with nonzero volume. For $`\lambda >0`$, define $`S_\mathrm{\Omega }(\lambda )`$ to be the sum of the values of $`T`$ at the isolated points of $`V`$ that belong to $`\lambda \mathrm{\Omega }`$. The mean value $`M_\mathrm{\Omega }`$ of $`T`$ over the isolated points of $`V`$ is defined by

$$M_\mathrm{\Omega }=\underset{\lambda \mathrm{}}{lim}\frac{S_\mathrm{\Omega }(\lambda )}{\mathrm{Vol}(\lambda \mathrm{\Omega })}.$$

### 4.2. Averaging a quasiperiodic trigonometric polynomial over the isolated points of a quasiperiodic semitrigonometric set

Consider the set of all frequencies of $`T`$ and $`T_{ij}`$. This set generates a subgroup $`𝒜`$ of $`(^n,+)`$, which is a direct sum of infinite cyclic subgroups:

$$𝒜=(𝒜_\mathcal{1})\mathrm{}(𝒜_𝒩).$$

That is, there exist $`A_1,\mathrm{},A_N^n`$ with no integral relations, such that each of the frequencies is an integral combination of $`A_1,\mathrm{},A_N`$.

Define a linear map $`\mathrm{\Phi }:^n^N`$ by $`\mathrm{\Phi }(x)=(A_1x,\mathrm{},A_Nx)`$. Then $`^n`$ acts on $`𝕋^N=^N/^N`$ by

$$x:\phi (\phi +\mathrm{\Phi }(x))\mathrm{mod}^N.$$

Let $`𝒪`$ be the orbit of this action through the origin. Since each of the frequencies $`\alpha _k`$ is an integral combination of $`A_1,\mathrm{},A_N`$, the functions $`T`$ and $`T_{ij}`$ are the restrictions to $`𝒪`$ of some trigonometric polynomials $`\stackrel{~}{T}`$ and $`\stackrel{~}{T}_{ij}`$ defined on $`𝕋^N`$. For example, $`\mathrm{cos}(2\pi \alpha x)`$ is the restriction of $`\mathrm{cos}2\pi m_i\phi _i`$, where $`\alpha =m_iA_i`$ for $`m_i`$.

To apply Theorem 3.4 here we need the orbit to be dense (Lemma 4.1), the mapping $`\mathrm{\Phi }`$ to be injective (Lemma 4.2), and the dimension $`n`$ of the orbit to be less than $`N`$ (if $`\mathrm{\Phi }`$ is injective then $`nN`$; $`N=n`$ corresponds to the periodic case which we treat in Remark 4.4).

###### Lemma 4.1.

There are no nonzero integral vectors orthogonal to the plane of the orbit $`𝒪`$. Therefore, $`𝒪`$ is dense in the torus $`𝕋^N`$.

###### Proof.

Let $`A`$ be the $`N\times n`$ matrix whose rows are the generators $`A_1,\mathrm{},A_N`$ of $`𝒜`$. We denote the columns of this matrix by $`A^1,\mathrm{},A^n`$. These columns generate the $`n`$-plane $`\mathrm{\Phi }(^n)`$ of the orbit $`𝒪`$. If there exists a nonzero integral vector $`(k_1,\mathrm{},k_N)`$ orthogonal to $`𝒪`$ then it is orthogonal to the vectors $`A^1,\mathrm{},A^n`$. This implies that $`k_1A_1+\mathrm{}+k_NA_N=0`$, i.e. there exists a nontrivial linear combination with integral coefficients on the generators $`A_1,\mathrm{},A_N`$ of $`𝒜`$, which gives a contradiction. By Appendix A the orbit $`𝒪`$ is dense in the torus $`𝕋^N`$. ∎

###### Lemma 4.2.

If the dimension of the orbit $`𝒪`$ is less than $`n`$, then the set $`V`$ described by the $`T_{ij}`$ has no isolated points, and the mean value $`M_\mathrm{\Omega }`$ is equal to zero.

###### Proof.

If the dimension of the orbit is less than $`n`$, then $`A^1,\mathrm{},A^n`$ are linearly dependent, and the solution space of the system $`A^1x_1+\mathrm{}+A^nx_n=0`$ is nontrivial. If $`T_{ij}(x)=0`$, then $`T_{ij}(x+y)=0`$ for each $`y`$ in the solution space. This means that $`V`$ has no isolated points, and the mean value is equal to zero. ∎

If the orbit is $`n`$dimensional and $`n<N`$, then $`^n`$ is mapped bijectively to a dense $`n`$-dimensional orbit on the real torus $`𝕋^N`$ and trigonometric polynomials $`T_{ij}`$ regarded as functions $`\stackrel{~}{T}_{ij}`$ on $`𝕋^N`$ define a set $`\stackrel{~}{V}`$. The intersection of the set $`\stackrel{~}{V}`$ and the orbit $`𝒪`$ is the semitrigonometric set $`V`$. We are adding up the values of $`\stackrel{~}{T}`$ over the isolated (in the topology of the orbit) points of intersection of the set $`\stackrel{~}{V}`$ and the orbit. Here is the main result of this chapter which now follows directly from Theorem 3.4.

###### Theorem 4.3.

Let $`V`$ be a quasiperiodic semitrigonometric subset of $`^n`$, $`T`$ a quasiperiodic trigonometric polynomial. Let $`\stackrel{~}{V}`$, $`\stackrel{~}{T}`$, $`𝒪`$ be the corresponding set, function on $`𝕋^N`$, and orbit, as constructed above. Let $`\mathrm{\Omega }`$ be a bounded subset of $`^n`$ with nonzero volume, whose boundary has fractal dimension less than $`n`$.

Then the mean value $`M_\mathrm{\Omega }`$ always exists and is equal to the $`𝒪`$-transversal$`\stackrel{~}{T}`$-weighted volume of $`\stackrel{~}{V}`$, assuming that the dimension of the orbit is $`n`$, and $`n<N`$. If the dimension of the orbit is less than $`n`$, the mean value is equal to zero.

###### Remark 4.4.

If the orbit has the same dimension as the torus, i.e., $`N=n`$, the set $`V^n`$ and the function $`T`$ have $`n`$ linearly independent periods $`A^1,\mathrm{},A^n`$, which define a torus $`^n/(A^1\mathrm{}A^n)`$. The mean value $`M_\mathrm{\Omega }`$ in this case is equal to the sum of the values of $`T`$ at the isolated points of $`V`$ in this torus, divided by the volume of the torus.

###### Remark 4.5.

This theorem holds and our proof works if $`V`$ is a quasiperiodic subanalytic subset of $`^n`$, that is, V is the intersection of a subanalytic subset $`\stackrel{~}{V}`$ of $`𝕋^N`$ with some dense $`n`$-dimensional orbit, and $`T`$ is the restriction to this orbit of some bounded subanalytic function on $`\stackrel{~}{V}`$.

###### Corollary 4.6.

If all the multiples $`k\alpha `$, $`k\{0\}`$ of each frequency $`\alpha `$ of $`T`$ do not belong to the subgroup $`𝒜`$ of $`(^n,+)`$ generated by the frequencies of $`T_{ij}`$, then the mean value $`M_\mathrm{\Omega }`$ is equal to zero.

###### Proof.

It is enough to check this statement for $`T`$ with a single frequency $`\alpha `$, i.e., for $`T(x)=c\mathrm{cos}2\pi \alpha x+d\mathrm{sin}2\pi \alpha x`$. The set of all frequencies of the trigonometric polynomials $`T_{ij}`$ describing the set $`V`$ generates a group $`𝒜^{}`$ in $`(^n,+)`$, which is the direct sum of infinite cyclic subgroups:

$$𝒜^{}=(𝒜_\mathcal{1})\mathrm{}(𝒜_{𝒩\mathcal{1}}).$$

Since no multiple of the frequency $`\alpha `$ of $`T`$ belongs to $`𝒜^{}`$, the group $`𝒜`$ obtained from $`𝒜^{}`$ by throwing in $`\alpha `$, is the direct sum:

$$𝒜=(𝒜_\mathcal{1})\mathrm{}(𝒜_{𝒩\mathcal{1}})(\alpha ).$$

Therefore, raising the situation to the torus $`^N/^N`$, we obtain a function $`\stackrel{~}{T}(\phi _N)=c\mathrm{cos}\phi _N+d\mathrm{sin}\phi _N`$ that depends only on the variable $`\phi _N`$, and a set $`\stackrel{~}{V}`$, which is described by the functions $`\stackrel{~}{T}_{ij}`$ that do not depend on $`\phi _N`$.

By Theorem 4.3, the mean value is equal to the $`𝒪`$-transversal $`\stackrel{~}{T}`$-weighted volume of $`\stackrel{~}{V}`$. Since the $`\stackrel{~}{T}_{ij}`$ do not depend on $`\phi _N`$, the manifold $`M_{Nn}(\stackrel{~}{V})`$ (which is the smooth $`Nn`$ dimensional, transversal to the orbit $`𝒪`$, part of $`\stackrel{~}{V}`$) is the direct product $`\stackrel{~}{V}^{}\times S^1`$, where $`\stackrel{~}{V}^{}`$ is a manifold in $`^{N1}/^{N1}`$, and $`S^1=/`$. The form $`\omega =p^{}d\phi _1\mathrm{}d\phi _N(\mathrm{},A^1,\mathrm{},A^n)`$ can be written as $`\omega _1d\phi _N`$ where $`\omega _1`$ is a form on $`\stackrel{~}{V}^{}`$. Therefore, the mean value is equal to

$$_{M_{Nn}(\stackrel{~}{V})}T\omega =_{\stackrel{~}{V}^{}}\left(_{S^1}(c\mathrm{cos}2\pi \phi _N+d\mathrm{sin}2\pi \phi _N)𝑑\phi _N\right)\omega _1=0.$$

## 5. Zeros of systems of exponential sums

In this section we apply our main result (Theorem 3.4) to the complex case, that is, to computing the mean value of an exponential sum over the zeros of a system of $`n`$ exponential sums in $`^n`$. We show that this mean value always exists. From the integral representation for the mean value that we obtain, we deduce a proof of Conjecture 2.3 in the case when the frequencies of the exponential sum are not commensurate with the frequencies of the system.

### 5.1. Definitions

Fix a system

(5.1)
$$F_1(z)=\mathrm{}=F_n(z)=0,z^n$$

of exponential sums with a developed collection of Newton polyhedra. Let $``$ be a family of systems of $`n`$ exponential sums with the same collection of spectra as (5.1), and whose coefficients have the same absolute values as the corresponding coefficients in (5.1). The following theorem is a particular case of a result proved by Gelfond .

###### Theorem 5.1.

There exists $`R>0`$ such that all the zeros of the systems from the family $``$ belong to the strip $`S_R\times \mathrm{Im}^\mathrm{n}`$, where $`S_R\mathrm{Re}^\mathrm{n}`$ is a ball of radius $`R`$ centered at the origin. This implies that all the zeros of the systems from $``$ are isolated.

Let $`G`$ be an exponential sum with real frequencies, and $`\mathrm{\Omega }`$ a bounded subset of $`\text{Im}^n`$ with nonzero volume. For $`\lambda >0`$, let $`S_\mathrm{\Omega }(\lambda )`$ be the sum of the values of $`G`$ at the zeros of (5.1) (counting multiplicities) that belong to the strip $`^n\times \lambda \mathrm{\Omega }^n.`$ Let $`S_\mathrm{\Omega }(\lambda )^k`$ be the sum of the values of $`G`$ at the zeros (not counting multiplicities) of the system (5.1) of multiplicity at least $`k`$ that belong to the strip $`^n\times \lambda \mathrm{\Omega }^n`$. By Theorem 5.1, $`S_\mathrm{\Omega }(\lambda )`$ and $`S_\mathrm{\Omega }(\lambda )^k`$ are well–defined.

Define

$$M_\mathrm{\Omega }=\underset{\lambda \mathrm{}}{lim}\frac{S_\mathrm{\Omega }(\lambda )}{\mathrm{Vol}(\lambda \mathrm{\Omega })},M_\mathrm{\Omega }^k=\underset{\lambda \mathrm{}}{lim}\frac{S_\mathrm{\Omega }^k(\lambda )}{\mathrm{Vol}(\lambda \mathrm{\Omega })}.$$

Theorem 5.2 below provides an integral formula for computing the mean values $`M_\mathrm{\Omega }`$, $`M_\mathrm{\Omega }^k`$ similar to the formula in the real case (Theorem 4.3). In particular, Theorem 5.2 implies that the mean value $`M_\mathrm{\Omega }`$ always exists. Another consequence of this integral formula is that in the case when the frequencies of $`G`$ are not commensurate with the frequencies of the exponential polynomials from the system, the mean value $`M_\mathrm{\Omega }`$ is equal to zero (Corollary 5.3).

### 5.2. Averaging an exponential sum over the zeros of a system of exponential sums

As in the real case, we deduce an integral formula for the mean value from Theorem 3.4. The exponential sums $`F_j`$, $`G`$ are quasiperiodic along the imaginary subspace of $`^n`$. Since the system (5.1) has a developed collection of Newton polyhedra, the real parts of the zeros of this system belong to some ball $`S_R\text{Re }^n`$. Therefore, we can restrict our attention to this ball, and think of the $`F_j`$ and $`G`$ as of functions quasiperiodic along the imaginary subspace and periodic along the real subspace. This allows us to consider $`^n`$ as an orbit on some real torus and apply Theorem 3.4 to this situation.

As before, the set of all frequencies $`\alpha `$ of the exponential sums $`F_j`$, $`G`$ generates a subgroup $`𝒜`$ in $`(^n,+)`$, which is a direct sum of infinite cyclic subgroups:

$$𝒜=(\alpha _\mathcal{1})\mathrm{}(\alpha _𝒩).$$

In other words, there exist $`\alpha _1,\mathrm{},\alpha _N^n`$ with no integral relations, such that each of the frequencies $`\alpha 𝒜`$ is an integral combination of $`\alpha _1,\mathrm{},\alpha _N`$.

Define a linear map

$$\mathrm{\Phi }_1:\text{Im}^n^N\text{by}\mathrm{\Phi }_1(y)=(\alpha _1y,\mathrm{},\alpha _Ny),y\text{Im}^n.$$

According to Theorem 5.1 the real parts of the zeros of the system (5.1) belong to the open ball $`S_R\text{Re}^n`$ centered at the origin. Let $`\mathrm{\Phi }_0:\text{Re}^n^n`$ be a composition of a shift and a rescaling that maps $`S_R`$ inside the unit cube. We now think of $`F_j`$ and $`G`$ as of functions periodic along the real subspace.

Then $`\mathrm{\Phi }=\mathrm{\Phi }_0\times \mathrm{\Phi }_1`$ maps $`^n=\text{Re }^n\times \text{Im }^n`$ to $`^n\times ^N`$. The map $`\mathrm{\Phi }`$ defines an action of $`^n`$ on $`𝕋^{n+N}`$ by

$$\phi (\phi +\mathrm{\Phi }(x,y))\mathrm{mod}^{n+N},\phi 𝕋^{n+N},x\text{Re }^n,y\text{Im }^n.$$

Let $`𝒪`$ be the orbit of this action through the origin. Since each of the frequencies $`\alpha `$ is an integral combination of $`\alpha _1,\mathrm{},\alpha _n`$, the sums $`F_j`$, $`G`$ are the restrictions to $`𝒪`$ of some functions $`\stackrel{~}{F}_j`$, $`\stackrel{~}{G}`$ defined on $`𝕋^{n+N}`$. We note that these functions may only be non-analytic at the points $`\phi `$ where at least one of the coordinates $`\phi _1,\mathrm{},\phi _n`$ is equal to zero.

The set of all zeros of the system (5.1) of multiplicity at least $`k`$ is an analytic set defined by exponential equations $`F_j(z)=0`$ and a few more exponential equations whose frequencies belong to the subgroup $`𝒜`$ of $`(^n,+)`$ (see, for example, ).

Let $`V^k`$ be the set of all zeros of the system of multiplicity greater than or equal to $`k`$. Consider the system of exponential sums that defines $`V^k`$. Regarded as functions on $`𝕋^{n+N}`$, these sums define a subset $`\stackrel{~}{V}^k`$ of the torus $`𝕋^{n+N}`$. We now state the main result of this section.

###### Theorem 5.2.

Consider a system of exponential sums with real frequencies

(5.2)
$$F_1(z)=\mathrm{}=F_n(z)=0,z^n$$

with a developed collection of Newton polyhedra. Let $`G`$ be an exponential sum with real frequencies, and $`\mathrm{\Omega }`$ a bounded subset of $`^n`$ with nonzero volume whose boundary has fractal dimension less than $`n`$. Let $`\stackrel{~}{V}^k`$, $`\stackrel{~}{G}`$, $`𝒪`$ be the corresponding subsets of $`𝕋^{n+N},`$ function on $`𝕋^{n+N}`$, and orbit, as constructed above.

Then if $`n<N`$ the mean value $`M_\mathrm{\Omega }^k`$ is equal to the $`\stackrel{~}{G}`$-weighted $`𝒪`$-transversal volume of $`\stackrel{~}{V}^k`$. The mean value $`M_\mathrm{\Omega }`$ is equal to the sum $`M_\mathrm{\Omega }^1+\mathrm{}+M_\mathrm{\Omega }^k`$ for some $`k`$.

###### Proof.

We first show that the dimension of the orbit $`𝒪`$ is $`2n`$. It is enough to show that the dimension of the plane $`\mathrm{\Phi }_1(\text{Im }^n)`$ is $`n`$. This plane is generated by the columns of the matrix of the linear map $`\mathrm{\Phi }_1`$. The rows of this matrix are the generators $`\alpha _1,\mathrm{},\alpha _N`$ of the group $`𝒜`$. Thus the dimension of the plane $`\mathrm{\Phi }_1(\text{Im }^n)`$ is equal to the dimension of the linear space generated by the frequencies $`\alpha _1,\mathrm{},\alpha _N`$. Since the collection of the Newton polyhedra of the system (5.2) is developed, its frequencies generate $`^n`$. We conclude that the dimension of the orbit is $`n`$.

Similarly to the real case, it is easy to see that there are no integral vectors orthogonal to the orbit $`𝒪`$. Therefore, the orbit $`𝒪`$ is dense in the torus $`𝕋^{n+N}`$ (see Lemma 4.1).

To derive Theorem 5.2 from Theorem 3.4 we show that $`\stackrel{~}{V}^k`$ is an analytic subset of $`𝕋^{n+N}`$ and $`\stackrel{~}{G}`$ is analytic on $`\stackrel{~}{V}^k`$. Each of $`\stackrel{~}{V}^k`$ is defined by a system of equations analytic everywhere on $`𝕋^{n+N}`$ but possibly at the points where at least one of the coordinates $`\phi _1,\mathrm{},\phi _n`$ is equal to zero. Notice that $`\stackrel{~}{V}^k`$ is contained in the “cylinder” $`\mathrm{\Phi }_0(S_R)\times 𝕋^N`$. Indeed, the intersections of $`\stackrel{~}{V}`$ with the orbit through the origin belong to this cylinder since the zeros of the initial system lie in the strip $`S_R\times \text{Im }^n`$. Restricting the functions $`\stackrel{~}{F}_j`$ to a shifted orbit we obtain functions that differ from $`F_j`$ by a shift along the imaginary subspace. By Theorem 5.1 the zeros of this new system belong to the strip $`S_R\times \text{Im }^n`$, and we conclude that $`\stackrel{~}{V}^k\mathrm{\Phi }_1(S_R)\times 𝕋^N`$. Therefore, $`\stackrel{~}{V}^k`$ is locally defined by a system of analytic equations, and $`G`$ is analytic on $`\stackrel{~}{V}^k`$.

To show that $`M_\mathrm{\Omega }`$ can be computed by the formula $`M_\mathrm{\Omega }=M_\mathrm{\Omega }^1+\mathrm{}+M_\mathrm{\Omega }^k`$ for some $`k`$, we need to prove that the sets $`V^k`$ are empty, starting from some $`k`$. Assume that there exists a sequence $`z_k=x_k+iy_k`$, $`x_kS_R`$, such that $`z_k`$ belongs to $`V^k`$. Raising this sequence to $`S_R\times 𝕋^N`$, we obtain a sequence $`\stackrel{~}{z}_k`$ in $`S_R\times 𝕋^N`$ such that $`\stackrel{~}{z}_k`$ belongs to $`\stackrel{~}{V}^k`$. Since $`S_R\times 𝕋^N`$ is compact, there exists a limit point $`\stackrel{~}{z}_0S_R\times 𝕋^N`$. This point $`\stackrel{~}{z}_0`$ belongs to each $`\stackrel{~}{V}^k`$. Although this point does not necessarily belong to the orbit $`S_R\times \mathrm{\Phi }_1(\text{Im }^N/^N)`$, it belongs to some shifted orbit. The restrictions of the functions $`\stackrel{~}{F}_j`$ to this shifted orbit define an exponential system whose zeros are isolated and belong to $`S_R\times \text{Im }^n`$ (Theorem 5.1). Therefore, $`\stackrel{~}{z}_0`$ is a zero of the shifted system of infinite multiplicity, which is impossible since all the zeros of this system are isolated. We proved that all $`V^k`$ are empty starting from some $`k`$.

Now the theorem follows from Theorem 3.4. ∎

The proof of the following corollary is a repetition of the proof of Corollary 4.6.

###### Corollary 5.3.

If all the multiples $`k\alpha ,k\{0\}`$ of each frequency $`\alpha `$ of $`F`$ do not belong to the subgroup $`𝒜`$ of $`(^n,+)`$ generated by the frequencies of the $`F_j`$, then the mean value $`M_\mathrm{\Omega }`$ is equal to zero.

## 6. Proof of Theorem 3.4

###### Proof.

Recall that $`𝒪=(𝒶+\mathcal{\Phi }(^𝓃))\mathrm{mod}^𝒩`$, where $`a𝕋^N`$, and $`\mathrm{\Phi }:^n^N`$ is a linear injective map. We send the torus $`𝕋^N`$ to the cube $`I=[0,1)^N^N`$ by cutting $`𝕋^N`$ along $`\phi _j=1`$. Choose a new basis $`\{B^1,\mathrm{},B^{Nn},A^1,\mathrm{},A^n\}`$ in $`^N`$ whose last $`n`$ coordinate vectors are the images of the standard basis in $`^n`$ under $`\mathrm{\Phi }`$, $`A^j=\mathrm{\Phi }(e_j)`$. Note that the vectors $`A^1,\mathrm{},A^n`$ generate the plane of the orbit $`\mathrm{\Phi }(^n)`$. Let the remaining vectors $`B^1,\mathrm{},B^{Nn}`$ form a basis of the $`(Nn)`$-plane $`L`$ orthogonal to the orbit, so that the determinant of the change of coordinates matrix from the standard basis in $`^N`$ to this new basis is equal to one.

Let $`x=(x_1,\mathrm{},x_n)`$ be the coordinates along the orbit, and $`y=(y_1,\mathrm{},y_{Nn})`$ the coordinates in the orthogonal plane $`L`$. We choose a cube $`Q=Q_y\times Q_x`$ whose edges go along the new coordinate axes and $`IQ`$. We now regard $`\stackrel{~}{V}`$ as a subset of $`IQ`$ and $`\stackrel{~}{T}(\phi )`$ as a function $`\stackrel{~}{T}(y,x)`$ on $`\stackrel{~}{V}I`$.

The set $`\stackrel{~}{V}`$ is a subset of $`Q`$, subanalytic in $`P^N`$, and the graph of the function $`\stackrel{~}{T}`$ is a subset of $`Q\times `$ subanalytic in $`P^{N+1}`$. Therefore, by the cell decomposition theorem (Theorem 8.3), there exists a decomposition of $`^N`$ with the fixed order of coordinates $`(y,x)`$ into finitely many cells such that $`\stackrel{~}{V}`$ is a union of cells and $`\stackrel{~}{T}`$ is $`C^1`$ on each of these cells.

Since the cell decomposition respects the standard projection, each cell is either a part of the closure of another cell, or is disjoint from it. Only $`(i_1,\mathrm{},i_{Nn},0,\mathrm{},0)`$-cells that do not belong to the closure of a cell of higher dimension may have isolated intersections with the orbit. We will show later that cells of dimension less than $`Nn`$ do not contribute to the mean value. Therefore, we can restrict our attention to $`(1,\mathrm{},1,0,\mathrm{},0)`$-cells with $`Nn`$ ones and $`n`$ zeros that do not belong to a closure of a cell of higher dimension.

Let $`C`$ be one of such cells. Then

$$C=\{(y,x)|yD,x=H(y)\},$$

where $`H=(H_1,\mathrm{},H_n)`$ for some $`C^1`$ functions $`H_1,\mathrm{},H_n`$, and $`D`$ is a $`(1,\mathrm{},1)`$-cell in $`Q_y`$.

Now we define a counting function $`g`$ on the torus $`𝕋^N`$. Through each point $`aC`$ we draw an $`n`$-plane parallel to the orbit. If $`\phi `$ belongs to $`(a+\mathrm{\Phi }(B_r))/^N`$, where $`B_r^n`$ is the ball of radius $`r`$ centered at the origin, we define $`g(\phi )=\stackrel{~}{T}(a)/v_r`$, where $`v_r`$ is the volume of the ball $`B_r^n`$. (Here $`r`$ is chosen so small that the mapping $`x\mathrm{\Phi }(x)/^N`$ restricted to $`B_r`$ is one-to-one.) Otherwise we define $`g(\phi )=0`$.

In $`(y,x)`$-coordinates, $`g`$ is defined by

$$g(y,H(y)+x)=\{\begin{array}{cc}\stackrel{~}{T}(y,H(y))/v_r,\hfill & \text{if }yD,xB_r\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$

Figure 2 illustrates the definition of $`g`$ in the case when the cell $`C`$ is one-dimensional. Through each $`aC`$ draw a vertical segment of length $`2r`$ centered at $`a`$. This segment is the ball $`a+B_r`$ from the general construction. Define $`g`$ to be identically equal to $`\stackrel{~}{T}(a)/2r`$ on that segment.

Since the function $`\stackrel{~}{T}`$ is continuous and bounded on $`C`$ and the functions defining $`C`$ are continuously differentiable, $`g`$ is Riemann integrable on the torus $`𝕋^N`$. We apply Weyl’s equidistribution law for multidimensional trajectories (Theorem 7.1):

(6.1)
$$\underset{\lambda \mathrm{}}{lim}\frac{1}{\mathrm{Vol}(\lambda \mathrm{\Omega })}_{\lambda \mathrm{\Omega }}g(\mathrm{\Phi }(x))𝑑x=_{𝕋^N}g(\phi )𝑑\phi .$$

Let us show first that the left hand side part of this identity is exactly the contribution of the cell $`C`$ to the mean value $`M_\mathrm{\Omega }`$. If $`a+B_r`$, a ball centered at the point $`a^n`$, lies entirely in $`\lambda \mathrm{\Omega }`$ then by the definition of $`g`$

$$_{a+B_r}g(\mathrm{\Phi }(x))𝑑x=T(a).$$

It remains to take care of the points $`a`$ such that $`a+B_r`$ does not lie entirely in $`\lambda \mathrm{\Omega }`$. The dimension of the boundary of $`\mathrm{\Omega }`$ is less than or equal to $`n\delta `$ for some $`\delta >0`$. Hence the number of balls of radius $`r`$ that are required to cover the $`r`$-neighborhood of the boundary of $`\lambda \mathrm{\Omega }`$ (i.e. the set of all points whose distance to the boundary is less than $`r`$) is bounded by $`c\lambda ^{n\delta }`$ for some constant $`c`$. The number of isolated points of $`V`$ in each of such balls is bounded by the number of cells in our cell decomposition. Therefore, the result of adding up $`T`$ over the isolated points of $`V`$, whose distance to the boundary of $`\lambda \mathrm{\Omega }`$ is less than $`r`$, is bounded by $`c_0\lambda ^{n\delta }`$ for some constant $`c_0`$. When we divide the result by the volume of $`\lambda \mathrm{\Omega }`$, which is proportional to $`\lambda ^n`$, the quotient approaches 0 as $`\lambda `$ approaches infinity. Therefore, we can disregard the points that are close to the boundary of $`\lambda \mathrm{\Omega }`$ while computing the mean value $`M_\mathrm{\Omega }`$. We have proved that

$$\underset{\lambda \mathrm{}}{lim}\frac{1}{\mathrm{Vol}(\lambda \mathrm{\Omega })}_{\lambda \mathrm{\Omega }}g(\mathrm{\Phi }(x))𝑑x$$

is the contribution of the cell $`C`$ to the mean value $`M_\mathrm{\Omega }`$.

We now work on the rewriting of the right hand side part of (6.1).

$`{\displaystyle _{𝕋^N}}g(\phi )𝑑\phi `$ $`=`$ $`{\displaystyle _Q}g(y,x)𝑑y𝑑x={\displaystyle _{Q_y}}\left({\displaystyle _{Q_x}}g(y,x)𝑑x\right)𝑑y`$
$`=`$ $`{\displaystyle _D}\stackrel{~}{T}(y,H(y))𝑑y={\displaystyle _C}\stackrel{~}{T}(\phi )p^{}𝑑\phi (\mathrm{},A^1,\mathrm{},A^n)={\displaystyle _C}\stackrel{~}{T}(\phi )p^{}\omega .`$
Next we note that cells of dimension less than $`Nn`$ do not contribute to the mean value. Indeed, for each of such cells we construct a counting function $`g`$ and repeat the above argument. It follows that the contribution of each of such cells is

$$_{𝕋^N}g(\phi )𝑑\phi ,$$

which is equal to zero since $`g`$ is zero on a set of full measure in $`𝕋^N`$.

The union of all $`(i_1,\mathrm{},i_{Nn},0,\mathrm{},0)`$-cells that do not belong to the closure of a cell of higher dimension is the $`(Nn)`$-dimensional manifold $`M_{Nn}(\stackrel{~}{V})`$, up to possibly a few cells of smaller dimension. Those cells of smaller dimension do not contribute to the mean value. We conclude that the mean value is equal to the $`\stackrel{~}{T}`$-weighted $`𝒪`$-transversal volume of $`\stackrel{~}{V}`$. ∎

## 7. Appendix A: Weyl’s equidistribution law for dense multidimensional orbits on the real torus.

Classical Weyl’s equidistribution law  states that the time-average of a Riemann integrable function along a one-dimensional dense orbit on a real torus coincides with the space average. Here we state a version of Weyl’s law for multidimensional orbits.

Let $`x^n`$, $`\phi 𝕋^N=^N/^N`$. A linear map $`\mathrm{\Phi }:^n^N`$ defines an action of $`^n`$ on $`𝕋^N`$ by $`\phi (\phi +\mathrm{\Phi }(x))\mathrm{mod}^N`$. Let $`𝒪`$ be an orbit of this action.

###### Theorem 7.1.

Let $`f`$ be a Riemann integrable function on the torus $`𝕋^N`$, $`𝒪`$ an orbit on $`𝕋^N`$ through an arbitrary point $`\phi _0𝕋^N`$ such that there are no nonzero integral vectors orthogonal to $`𝒪`$, and $`\mathrm{\Omega }`$ a subset of $`^n`$ with nonzero volume.

Then

$$\frac{1}{\mathrm{Vol}(\lambda \mathrm{\Omega })}_{\lambda \mathrm{\Omega }}f(\phi _0+\mathrm{\Phi }(x))𝑑x_{𝕋^N}f(\phi )𝑑\phi ,\text{as}\lambda \mathrm{}.$$

The proof of this theorem in  is a direct generalization of Weyl’s original argument for one-dimensional orbits.

###### Remark 7.2.

This theorem implies that an orbit is dense if and only if there are no nonzero integral vectors orthogonal to the plane of the orbit.

## 8. Appendix B: Cell-decomposition theorem

The cell-decomposition theorem for subanalytic sets states that a subanalytic set can be partitioned into finitely many cells, which are subanalytic subsets of especially simple form. We will follow the formulation of the cell-decomposition theorem as presented in .

###### Definition 8.1.

Let $`(i_1,\mathrm{},i_n)`$ be a sequence of zeros and ones of length $`n`$. An $`(i_1,\mathrm{},i_n)`$-cell is a subset of $`^n`$ obtained by induction on $`n`$ as follows:

* a (0)-cell is a point in $``$, a (1)-cell is an interval $`(a,b)`$;
* suppose $`(i_1,\mathrm{},i_n)`$-cells are already defined; then an $`(i_1,\mathrm{},i_n,0)`$-cell is the graph $`\mathrm{\Gamma }(f)`$ of a $`C^1`$ function $`f`$ on an $`(i_1,\mathrm{},i_n)`$-cell provided that $`\mathrm{\Gamma }(f)^{n+1}`$ is subanalytic in $`P^{n+1}`$; further, an $`(i_1,\mathrm{},i_n,1)`$-cell is a set

$$(f,g)_X:=\{(x,r)X\times :f(x)<r<g(x)\},$$

where $`X`$ is an $`(i_1,\mathrm{},i_n)`$-cell, functions $`f,g`$ are $`C^1`$ functions on $`X`$ whose graphs $`\mathrm{\Gamma }(f),\mathrm{\Gamma }(g)^{n+1}`$ are subanalytic in $`P^{n+1}`$, and $`f<g`$ on $`X`$. The constant functions $`f(x)=\mathrm{}`$, $`g(x)=+\mathrm{}`$ are also allowed.

###### Definition 8.2.

A decomposition of $`^n`$ is a special kind of partition of $`^n`$ into finitely many cells. The definition is by induction on $`n`$:

* a decomposition of $``$ is a collection

$$\{(\mathrm{},a_1),(a_1,a_2),\mathrm{},(a_k,+\mathrm{}),\{a_1\},\mathrm{},\{a_k\}\}$$

where $`a_1<\mathrm{}<a_k`$ are points in $``$.
* a decomposition of $`^{n+1}`$ is a finite partition into cells such that the set of projections of these cells to the first $`n`$ coordinates forms a decomposition of $`^n`$.

###### Theorem 8.3.

(The cell decomposition theorem) Let $`X`$ be a subset of $`^n`$ subanalytic in $`P^n`$. Let $`f`$ be a function on $`X`$ whose graph in $`^{n+1}`$ is a subanalytic subset of $`P^{n+1}`$. Then there exists a (finite) decomposition $`𝒟`$ of $`^n`$ such that X is a union of cells in $`𝒟`$ and the restriction $`f|_B:B`$ to each cell $`B𝒟`$ is $`C^1`$.


format: text

text:
# Structural, orbital, and magnetic order in vanadium spinels

## Abstract

Vanadium spinels (ZnV2O4, MgV2O4, and CdV2O4) exhibit a sequence of structural and magnetic phase transitions, reflecting the interplay of lattice, orbital, and spin degrees of freedom. We offer a theoretical model taking into account the relativistic spin-orbit interaction, collective Jahn-Teller effect, and spin frustration. Below the structural transition, vanadium ions exhibit ferroorbital order and the magnet is best viewed as two sets of antiferromagnetic chains with a single-ion Ising anisotropy. Magnetic order, parametrized by two Ising variables, appears at a tetracritical point.

Frustrated magnetic systems have generated interest among researchers for more than a decade Schiffer96 ; Greedan01 ; Moessner01 . Frustration refers to the inability of spins to satisfy conflicting interactions resulting in a large degeneracy of the ground state at the classical level. The effect is particularly strong for spins with antiferromagnetic interactions on the three-dimensional network of corner-sharing tetrahedra, known as the pyrochlore lattice Moessner98 . Magnets of this kind have rather unusual properties. For example, spinels containing magnetic ions Cr3+ (ZnCr2O4, MgCr2O4, and CdCr2O4) undergo a magnetoelastic phase transition Lee00 ; Tch02 , even though they have no spin chains commonly associated with the spin-Peierls effect Bray75 .

In this Letter I discuss magnetic and structural properties of another remarkable group of magnetic spinels containing V3+ ions (ZnV2O4, MgV2O4, and CdV2O4). These compounds represent a step up in complexity from their chromium counterparts because vanadium ions have an orbital degree of freedom—in addition to lattice vibrations and spin. The interplay of orbital, spin, and vibrational motion against the backdrop of a high spatial symmetry leads to a sequence of phase transitions involving changes in structural and magnetic properties. The physics involves the relativistic spin-orbit interaction, collective Jahn-Teller effect, and spin frustration.

The vanadium spinels undergo a structural phase transition at a temperature $`T_s`$ and a magnetic ordering at a lower temperature $`T_N`$. The structural distortion lowers the crystal symmetry from cubic $`Fd\overline{3}m`$ to tetragonal $`I4_1/amd`$ Onoda03 ; Reehuis03 , with a slightly shorter lattice period along the $`c`$ axis. The Néel temperature $`T_N`$ is rather small compared to the Curie-Weiss scale $`\mathrm{\Theta }_{\mathrm{CW}}`$, a sure sign of strong spin frustration. Ordered magnetic moments point along the $`c`$ axis and have the length $`\mu <\mu _B`$ at helium temperatures. Magnetic susceptibility is history-dependent both below and above $`T_N`$ Ueda97 ; Reehuis03 . The numbers for ZnV2O4 are representative of the entire family: $`c/a=0.994`$, $`T_s=51`$ K Ueda97 , $`T_N=42`$ K Mamiya97 , $`\mathrm{\Theta }_{\mathrm{CW}}=420`$ K Mamiya95 , and $`\mu =0.63\mu _B`$ Reehuis03 .

The main difference between Cr3+ and V3+ ions is the number of electrons in the partially filled $`3d`$ shell. The crystal field splits the $`3d`$ quintet into a low-energy triplet $`t_{2g}`$ and a high-energy doublet $`e_g`$, separated by a gap of roughly 2 eV. In accordance with Hund’s rules, the three $`d`$ electrons of Cr3+ have spin $`S=3/2`$ and occupy all three triplet levels. The ion has a spin degree of freedom only. In contrast, the two $`d`$ electrons of V3+ have spin $`S=1`$ and occupy two out of three $`t_{2g}`$ orbitals.

Tsunetsugu and Motome Tsunetsugu03  have recently developed a model explaining structural and magnetic properties of the vanadium spinels. Following the approach of Kugel and Khomskii, they focus on Coulomb and exchange interactions between magnetic ions. However, their model is at odds with experimental observations: the orbital order predicted by the theory is incompatible with the spatial symmetry of the tetragonal phase, identified as $`I4_1/amd`$ Onoda03 ; Reehuis03 . The staggering of ions with empty $`yz`$ and $`zx`$ orbitals would break reflections $`m`$ in the planes $`(110)`$ and $`(1\overline{1}0)`$ and diamond glides $`d`$ in the planes $`(100)`$ and $`(010)`$, lowering the symmetry to $`I4_1/a`$ \[Fig. 1(a)\].

I propose an alternative theory that is compatible with the structural data and predicts rather nontrivial magnetic behavior. In the new model, the structural order arises from simple single-ion physics: the interplay of the JT effect and relativistic spin-orbit interaction. The $`t_{2g}`$ orbitals are partially ordered: their occupation numbers satisfy $`n_{yz}=n_{zx}<n_{xy}`$ \[Fig. 1(b)\]. The distorted magnet is best decribed as a collection of antiferromagnetic chains, running in the directions $`[110]`$ and $`[1\overline{1}0]`$, with an effective spin $`J=2`$ and a single-ion anisotropy of the Ising type. Parallel chains are weakly coupled; the coupling of orthogonal chains is frustrated by symmetry. As a result, the Néel phase is characterized by two independent Ising order parameters; the magnetic transition occurs at a tetracritical point. Two nearby magnetic phases exist in which only one half of the chains order. Finally, memory effects in magnetic susceptibility are tied to the spontaneous breaking of the lattice symmetry.

Choice of the model. The main difficulty in obtaining a consistent theoretical picture of the vanadium spinels is a poor separation of energy scales. The $`t_{2g}`$ orbitals of a B site in spinels point away from the oxygen ligands and thus have a relatively weak JT coupling Goodenough63  comparable in magnitude to the relativistic spin-orbit coupling. (A recent estimate Mizokawa96  puts the JT energy in YVO3 at 20–60 meV. The spin-orbit splitting of the $`{}_{}{}^{3}F_{2}^{}`$ and $`{}_{}{}^{3}F_{3}^{}`$ levels in a free V3+ equals 39 meV Moore71 .) Furthermore, the effective interaction of spin and orbital degrees of freedom on adjacent vanadium ions has a similar energy scale Tsunetsugu03 . It is not at all obvious which of the three interactions plays the dominant role.

To make progress, one can adopt a phenomenological approach and investigate three simple limits, in which the dominant interactions are, respectively, (I) the JT coupling, (II) the spin-orbit coupling, (III) the Kugel-Khomskii interaction. Model III Tsunetsugu03  yields orbital order incompatible with the observed symmetry of the tetragonal phase Onoda03 ; Reehuis03 . We therefore explore limits I and II.

Basis states and notation. If the energy difference between the $`{}_{}{}^{3}P`$ and $`{}_{}{}^{3}F`$ electronic terms of the $`3d^2`$ configuration ($`1`$ eV) is neglected in comparison with its crystal-field splitting in a VO6 octahedron ($`2`$ eV), the $`{}_{}{}^{3}T_{1g}^{}`$ ground state of a V3+ ion has the $`t_{2g}^2e_g^0`$ configuration. Its orbital part is triply degenerate and has the following basis:

$`|X`$ $`=`$ $`(|zx|xy|xy|zx)/\sqrt{2},`$
$`|Y`$ $`=`$ $`(|xy|yz|yz|xy)/\sqrt{2},`$ (1)
$`|Z`$ $`=`$ $`(|yz|zx|zx|yz)/\sqrt{2}.`$
In this Hilbert space, the operator of orbital momentum $`𝐋`$ can be expressed in terms of an effective angular momentum $`𝐋^{}`$ of length $`L^{}=1`$ such that

$$L_z^{}|X=i|Y,L_z^{}|Y=i|X,L_z^{}|Z=0$$
(2)
etc. It can be verified that $`𝐋=\alpha 𝐋^{}`$, where $`\alpha =1`$ in the approximation made (1). Inclusion of electronic correlations mixes in the states $`t_{2g}^1e_g^1`$ and $`t_{2g}^0e_g^2`$ and puts $`\alpha `$ somewhere in the range between 1 and 3/2 Abragam51 .

Model I. A tetragonal distortion of a VO6 octahedron along the $`c`$ axis splits the $`{}_{}{}^{3}T_{1g}^{}`$ state into a doublet and a singlet. The splitting can be mimicked by the traceless operator

$$V_{\mathrm{JT}}=\delta [L_{}^{}{}_{z}{}^{2}L^{}(L^{}+1)/3],$$
(3)
where $`\delta `$ is the tetragonal strain expressed in energy units. For a $`3d^2`$ ion, $`\delta >0`$ for an elongated octahedron Goodenough63 . The lowest-energy state is $`|Z`$, in which the $`xy`$ orbital is empty and the energy is $`2\delta /3`$. In the opposite case (a flattened octahedron, $`\delta <0`$), the orbital ground state is doubly degenerate and has a higher energy $`|\delta |/3`$. Inclusion of the elastic energy $`𝒪(\delta ^2)`$ and minimization with respect to $`\delta `$ yields a ground state with $`\delta >0`$. Thus Model I, dominated by the JT coupling, predicts an elongation of VO6 tetrahedra contrary to the data Onoda03 ; Reehuis03 ; Ueda97 ; Mamiya97 ; Lee03 .

Model II. The relativistic spin-orbit coupling in V3+

$$V_{LS}=\lambda (𝐋𝐒)=\alpha \lambda (𝐋^{}𝐒)$$
(4)
has the strength $`\lambda +20`$ meV Mizokawa96 . The “total angular momentum” $`𝐉^{}=𝐋^{}+𝐒`$ is a conserved quantity. The lowest-energy state is the $`J^{}=2`$ quintuplet separated by the energy gap $`2\alpha \lambda `$ from the $`J^{}=1`$ levels.

Inclusion of the JT interaction (3) as a perturbation ($`\delta \lambda `$) splits the $`J^{}=2`$ quintuplet into three levels. With the aid of the Wigner-Eckart theorem, one finds

$$V_{\mathrm{JT}}=(\delta /2)[J_{}^{}{}_{z}{}^{2}J^{}(J^{}+1)/3].$$
(5)
Remarkably, the ground state has the same energy $`|\delta |`$ regardless of the sign of the $`\delta `$. In an elongated VO6 octahedron ($`\delta >0`$), the ground state is nondegenerate and has $`J_z=0`$; its magnetic moment $`\mu _B(\alpha 𝐋^{}+2𝐒)`$ is quenched. In a flattened octahedron ($`\delta <0`$), the ground state is a non-Kramers doublet $`J_z^{}=L_z^{}+S_z=\pm 2`$. Inclusion of the elastic energy $`𝒪(\delta ^2)`$ fails to determine the sign of the tetragonal JT distortion. It will be determined by several competing perturbations, such as higher-order JT terms, Coulomb and spin-exchange interactions between vanadium ions.

The effect of Coulomb interaction can be crudely estimated by counting the number of filled $`t_{2g}`$ orbitals with the largest overlap. For example, the Coulomb energy of a  V–V bond is raised if the $`xy`$ orbitals are occupied in both ions Tsunetsugu03 . In a state with uniform orbital order, the average occupations of the orbitals are as follows:

$$\begin{array}{ccc}J_z^{}=\pm 2:\hfill & n_{yz}=n_{zx}=1/2,\hfill & n_{xy}=1,\hfill \\ J_z^{}=0:\hfill & n_{yz}=n_{zx}=5/6,\hfill & n_{xy}=1/3.\hfill \end{array}$$
(6)
The expectation value of the Coulomb energy, summed over nearest-neighbor V–V pairs, will be proportional to $`n_{yz}^2+n_{zx}^2+n_{xy}^2`$. This expression evaluates to 3/2 for both $`J_z^{}=0`$ and $`J_z^{}=\pm 2`$ ground states. Thus the Coulomb repulsion between $`V^{3+}`$ ions is the same in the elongated and flattened ground states of VO6 octahedra and does not influence the outcome of their competition.

(The accidental degeneracy of the Coulomb term can be traced to the degeneracy of the flattened and elongated JT states for $`J^{}=2`$. Both the Coulomb and JT energies can be expressed in terms of the operators $`J_{}^{}{}_{z}{}^{2}J^{}(J^{}+1)/3`$ and $`J_{}^{}{}_{x}{}^{2}J_{}^{}{}_{y}{}^{2}`$. It is therefore likely that the conclusion reached in the previous paragraph is not sensitive to the approximations made in this paper.)

Exchange energy, on the other hand, clearly favors the magnetic ground states $`J_z^{}=\pm 2`$ because it can be made negative by an appropriate choice of spin orientations. In contrast, exchange energy is zero if the ions are in the nonmagnetic state $`J_z^{}=0`$. Anharmonic terms of the JT energy $`𝒪(\delta ^3)`$ can favor either sign of $`\delta `$. Therefore, one cannot deduce the sign of the distortion on general theoretical grounds. Nonetheless, experimental observation of antiferromagnetic order in ZnV2O4 at low temperatures points to the magnetically active states $`J_z^{}=\pm 2`$ characterized by a flattening of the VO6 octahedra along the $`c`$ axis. For a given direction of the spin, $`S_z=\pm 1`$, one of the electrons occupies the $`xy`$ orbital, whereas the other is in the state $`|yz\pm i|zx`$, whose actual orbital moment $`L_z=\alpha L_z^{}=\alpha `$ points opposite to the spin. Equal occupation of the $`yz`$ and $`zx`$ orbitals in such a state is in agreement with the observed spatial symmetry of the crystal ($`I4_1/amd`$). Magnetic moments of length $`(2\alpha )\mu _B`$, where $`1\alpha <3/2`$, point along the $`c`$ axis, as found experimentally. Thus phenomenological Model II, with strong relativistic spin-orbit coupling, a moderate JT distortion, and weak V–V interactions, is compatible with experimental observations.

Structural order and symmetry-breaking fields. Below $`T_s`$ the VO6 octahedra flatten along one of the three major axes. From the symmetry viewpoint, the long-range order is similar to that of the 3-state Potts ferromagnet. A two-component order parameter $`𝐟=(f_1,f_2)`$ can be constructed, e.g. from the lattice constants:

$$f_1=(a+b2c)/\sqrt{6},f_2=(ab)/\sqrt{2}.$$
(7a)
Other observables can be used to define the order parameter. The lattice constants $`(a,b,c)`$ can be replaced with orbital populations $`(n_{yz},n_{zx},n_{xy})`$, or with magnetic susceptibilities $`\chi _{ij}=M_j/H_i`$:
$$f_1=(\chi _{xx}+\chi _{yy}2\chi _{zz})/\sqrt{6},f_2=(\chi _{xx}\chi _{yy})/\sqrt{2}.$$
(7b)
Explicit construction of the order parameter allows us to identify symmetry-breaking fields that can be applied to select a particular phase with broken symmetry. The most obvious among these is the two-component (tetragonal and orthorombic) stress $`𝐡=(h_1,h_2)`$:

$$h_1=(e_{xx}+e_{yy}2e_{zz})/\sqrt{6},h_2=(e_{xx}e_{yy})/\sqrt{2}.$$
(8a)
Cooling under a weak uniaxial stress applied along the $`c`$ axis will produce an ordered phase with $`a=b>c`$.
A less obvious symmetry-breaking perturbation is magnetic field $`𝐇`$. The magnetic energy density of a sample can be written as $`(\chi _{xx}H_x^2+\chi _{yy}H_y^2+\chi _{zz}H_z^2)/2.`$ Unless the field is parallel to a $`111`$ direction, it lowers the energies of different ordered phases by different amounts. The symmetry-breaking field is

$$h_1=(H_x^2+H_y^22H_z^2)/\sqrt{6},h_2=(H_x^2H_y^2)/\sqrt{2}.$$
(8b)
These considerations could explain memory effects in magnetic susceptibility of ZnV2O4 above the Néel temperature $`T_N`$ Ueda97 ; Reehuis03 . A crystal cooled through the structural transition in magnetic field chooses the distorted phase with the lowest magnetic energy. The difference between $`\chi _{\mathrm{FC}}`$ and $`\chi _{\mathrm{ZFC}}`$ in powder samples may simply reflect the presence of the order parameter (7b) below $`T_s`$ and its sensitivity to the symmetry-breaking field (8b). It would be highly desirable to study these memory effects in single crystals.

Aniferromagnetic Ising chains with frustrated coupling. In the distorted state, the magnetic system, comprising spins and orbital moments, can be thought of as a collection of chains running in the $`[110]`$ and $`[1\overline{1}0]`$ directions (Fig. 2). Because the $`xy`$ orbitals are always occupied, the largest exchange coupling takes place along the chains. Exchange coupling along the other $`110`$ directions is expected to be substantially weaker because the relevant orbitals ($`yz`$ and $`zx`$) are only half-filled. Onoda and Hasegawa Onoda03  have extracted the values of $`J_{[110]}5`$ meV and $`J_{[011]}=J_{[101]}3`$ meV for CdV2O4. In addition to the single-ion Ising anisotropy (5), which is a product of orbital order and spin-orbit coupling, the exchange interactions are likely anisotropic as well.

As the distorted magnet is cooled further, spin correlations develop first along the chains, as indeed observed Lee03 . If the interchain coupling is weak, long-range magnetic order will set in at a low temperature, when each chain already has a large correlation length. In that case, the magnetic transition will have almost no effect on the energetics of the spins. Indeed, measurements of specific heat in MgV2O4 reveal a very weak singularity at the magnetic transition Mamiya97 .

While magnetic interactions between parallel chains are expected to be weak because they are well separated in space, the coupling between crossing chains may be substantial: four of their spins are nearest neighbors forming a tetrahedron (Fig. 2). However, in this case the coupling is frustrated by the lattice symmetry: a (110) reflection $`m`$ (or a glide $`a`$) changes the sign of staggered magnetization on the chain running in the direction $`[1\overline{1}0]`$ but preserves staggered magnetization on the  chain. Therefore, even if the crossing-chain coupling is numerically comparable to the intrachain coupling, it is still rendered weak by geometrical frustration.

Magnetic order. Neutron scattering indicates that adjacent parallel chains have equal values of staggered magnetization (Fig. 2). Because of the geometrical frustration, the system has four ground states: those shown in Fig. 2 and their time-reversed copies. These ground states can be characterized by two independent Ising order parameters: one staggered magnetization for each set of chains, $`\sigma _{[110]}`$ and $`\sigma _{[1\overline{1}0]}`$ chains-Tsu .

Magnetic tetracritical point. Both sets of chains appear to be ordered at the lowest temperatures Reehuis03 . In the simplest scenario, both Ising order parameters develop continuously at the magnetic transition temperature $`T_N`$, which is characteristic of a tetracritical point Fisher74 ; Rohrer77 . One should be able to observe two additional phases in the vicinity of the tetracritical point, in which only one set of chains is ordered, while the other is not. This may be achieved by creating an orthorombic strain $`e_{xy}`$ in the vicinity of the magnetic transition. Doing so will violate the equivalence of the  and $`[1\overline{1}0]`$ chains while still preserving the geometrical frustration. A temperature scan under a uniaxial  stress should reveal two separate magnetic transitions for each set of chains. The temperature–strain phase diagram in the vicinity of the tetracritical point is shown in Fig. 3.

Outlook. A phenomenological approach taken in this paper is justified by the complexity of vanadium spinels: the relevant degrees of freedom include lattice vibrations, orbital motion, and spin. It is desirable to put these considerations on a quantitative footing with the aid of a microscopic theory. Nonetheless, the qualitative predictions of this work—a spin gap due to single-ion anisotropy, a magnetic tetracritical point, and a potential relation between magnetic memory effects and broken lattice symmetry—can be tested experimentally.

Upon completion of this work I learned that a similar model had been considered by D.I. Khomskii.

Acknowledgment. I thank C. Broholm, D.I. Khomskii, A. Krimmel, M.Z. Hasan, S.-H. Lee, D. Louca, M. Onoda, P.G. Radaelli, O.A. Starykh, and H. Ueda for useful discussions. This work was supported in part by the U.S. National Science Foundation under the Grant No. DMR-0348679.


text:
# Numerical Investigation of the Entropy Crisis in Model Glass Formers

## I Introduction

The nature of the glass transition is one of the great unsolved problems in condensed matter physics. The kinetic glass transition (occurring at temperature $`T_g`$) is phenomenologically defined to occur when a liquid reaches a threshold viscosity Angell . It has long been speculated that a thermodynamic transition to a disordered solid would occur if the liquid could be cooled infinitely slowly without crystallizing. Kauzmann Kauzmann  argued that extrapolating the supercooled liquid entropy into the glassy regime implies that a temperature (known as the “Kauzmann temperature”, $`T_K`$) exists where the configurational entropy of the liquid becomes less than that of the crystal. Gibbs and coworkers Gibbs65  resolved this crisis by suggesting that a phase transition to an “ideal glass” occurs at $`T_K`$. Many theories of glassy behavior contain thermodynamic transitions to an ideal glass state. Gibbs and Di Marzio GibbsDimarzio  postulated that a second order transition occurs at $`T_K`$ while theories based upon mean-field spin glass ideas contain a random first order transition at $`T_K`$ XiaKirkKirk2b ; ParisiThermoParisiThermoLJ ; ParisiThermoSoft . Stillinger has argued that localized defects must destroy a second order ideal glass transition Stillinger . In general, localized excitations (non-mean-field effects) should round any transition that would occur at $`T_K`$, and recent arguments have been put forward to access the conditions under which the signatures of mean-field behavior survive Eastwood . While the above theories in some sense presuppose the existence of a transition at $`T_K`$ that drives the formation of a disordered solid, other theoretical models of glassy behavior completely avoid the notion of an underlying thermodynamic transition. In particular, a large class of kinetic facilitated models successfully explain many of the features exhibited by supercooled liquids without invoking the notion of configurational entropy Fredrickson1Fredrickson2Butler1JackleRitort . The difficulty of directly probing slowly cooled, low temperature glassy configurations, either experimentally or computationally, has made resolution of the issue of the entropy crisis extremely difficult.

Advances in simulation techniques have made the low temperature equilibration of certain model systems possible. These techniques include parallel tempering Marinari92Falcioni99Neirotti00 , histogram re-weighting Yan  and the use of specific non-local MC moves Gazzillo ; Krauth ; Swapping . Recent studies of model systems have produced seemingly conflicting results regarding the possible existence of an entropy crisis. Grigera and Parisi have demonstrated that a soft-sphere mixture displays a distinct signature of a thermodynamic transition near the Kauzmann temperature predicted by the random first order phase transition scenario Swapping , while Santen and Krauth have shown that a two dimensional hard-sphere system displays no such signature Krauth . It has been difficult to obtain further insight into these phenomena, as the simulation methodologies are highly system specific. Although extremely useful for several systems Krauth ; YukawaSpheres ; KrauthSanten ; MoreKrauth , the cluster algorithm of Krauth Cluster  has been hard to generalize to more complex and three dimensional systems, while the swapping method used by Grigera and Parisi fails for systems with attractive interactions. Other methods, such as simulated annealing and parallel tempering, have been found to equilibrate slowly in the glassy regimeMichele .

In this paper, we re-examine claims in favor and against the existence of an entropy crisis at a temperature $`T_K`$. We carefully explore the nature of the thermodynamic signature in the heat capacity of the soft-sphere system studied by Grigera and Parisi and study thermalized versions of system studied by Santen and Krauth in both two and three dimensions. A critical assessment is made concerning the ability of various non-local Monte Carlo techniques to equilibrate liquids under conditions necessary to probe the putative entropy crisis. Our conclusions suggest that strong claims for numerical evidence both in favor and against an entropy crisis in all systems studied to date must be re-evaluated.

## II Computational Details

When a liquid is deeply supercooled, its microscopic motion is characterized by rapidly growing relaxation times. This slow relaxation renders the use of conventional molecular dynamics and Monte Carlo techniques impractical for obtaining well equilibrated configurations near the glass transition temperature. Recently, several methodologies have been introduced in an effort to overcome the problems associated with the rapid onset of sluggish dynamics that hinder equilibration with local simulation techniques. Parallel tempering methods have been shown to speed up equilibration above the mode-coupling temperature $`T_c`$ YamaKob , but have not been useful in obtaining well equilibrated configurations significantly below this temperature Michele . Other techniques, such as expanded ensemble methods PabloYan , multi-canonical algorithms Berg , and the density-of-states Monte Carlo approach WangWang2 ; Yan2Pablo2  may provide routes to obtaining equilibrated configurations at low temperatures, although few systematic studies have been performed for glassy systems with these techniques and we will not discuss them further here Faller ; Note .

Conventional Monte Carlo techniques augmented with non-local moves offer a powerful approach to overcoming ergodicity problems associated with rough free energy landscapes. These methods allow for moves that interchange particles or groups of particles while still maintaining detailed balance. Even if such moves are rarely accepted, they greatly accelerate sampling by providing routes for relaxation (e.g. the hopping of a particle out of the cage comprised of its neighbors) that would take many local Monte Carlo moves to achieve. Recently, two types of non-local Monte Carlo methods have been applied to the study of glassy systems: The cluster pivot algorithm (CPA) of Krauth and coworkers Krauth  and the swap Monte Carlo (SMC) algorithm applied by Grigera and Parisi to the study of the entropy crisis Swapping .

In this work we introduce a trivial modification of the SMC approach, called swap-sector Monte Carlo (SSMC), that is optimized for equilibrating polydisperse systems. As with SMC, the positions of particles of different sizes are swapped at varying intervals. Noting that the probability of swapping falls off quickly as the ratio of particle sizes differs from 1, we introduce a swap-sector, where particles are only swapped with those that differ in radius by less than a tunable parameter $`\mathrm{\Delta }\sigma `$. In the limit of large $`\mathrm{\Delta }\sigma `$, standard SMC is recovered. $`\mathrm{\Delta }\sigma `$ can be optimized for efficiency, as a large $`\mathrm{\Delta }\sigma `$ yields a low acceptance ratio, and a small $`\mathrm{\Delta }\sigma `$ allows pairs of particles to be swapped back and forth repeatedly. This method has been used to equilibrate both hard and soft-sphere polydisperse systems, facilitating direct comparison between the two.

Three different systems are simulated in this work. The first is a binary soft-sphere system identical to that studied by Grigera and Parisi, with

$$V=\underset{ij}{\overset{N}{}}\left(\frac{\sigma _i+\sigma _j}{|𝐫_𝐢𝐫_𝐣|}\right)^{12},$$
(1)
where $`\sigma _1/\sigma _2=1.2`$, $`\sigma _{12}=\frac{1}{2}(\sigma _1+\sigma _2)`$ and $`ϵ_1=ϵ_2=1`$. The units are set by choosing $`\sigma _1=1`$. Runs are executed by varying the temperature with constant volume, chosen so that $`\rho ^{}=\rho \sigma _x^3=1`$, where

$$\sigma _x^3=\frac{1}{4}\sigma _1^3+\frac{1}{2}\sigma _{12}+\frac{1}{4}\sigma _2^3$$
(2)
is the effective radius defined by conformal solution theory HansenMcDonald . All simulations in this work employ periodic boundary conditions. Systems of 34, 60, 70 and 258 particles are studied, consistent with previous studies. Although clearly too small to study dynamics in a regime where the dynamical correlation length is noticeably growing, our studies should provide a useful probe of the underlying thermodynamic behavior.

The second system is the two dimensional hard-sphere system of Santen and Krauth Krauth , run in the NPT ensemble. The 256 particles have radii ranging from $`R_1=1`$ to $`R_{256}=20`$, with $`R_iR_{i1}=\mathrm{\Delta }`$ for constant $`\mathrm{\Delta }`$. Isothermal compressibilities are calculated both by finite difference of the volume with respect to pressure, $`\kappa _{diff}=\frac{1}{V}\frac{V}{P}`$, and by fluctuations of the volume, $`\kappa _{fluc}=\frac{V^2V^2}{V}`$. If the system is well equilibrated and the data well averaged, $`\kappa _{diff}=\kappa _{fluc}`$. This provides a necessary (but not sufficient) test for equilibration that, together with convergence and stability of thermodynamic averages, is applied throughout. We also applied the Flyvberg-Petersen method Flyvberg  to determine accurate estimates of the statistical error, which were well below 0.1% of the heat capacity, even at extremely low temperatures.

Lastly, a polydisperse soft-sphere system is studied with $`\sigma _i=1+\mathrm{\Delta }i`$ where $`1iN`$ and $`\mathrm{\Delta }`$ is chosen so that $`R_{max}=20,`$ providing a thermal version of the Santen and Krauth system. Runs are executed by both varying the temperature with constant volume, and varying the volume with constant temperature to confirm the validity of the results (excess thermodynamic properties in soft-sphere systems depend only on the parameter $`\mathrm{\Gamma }=\rho T^{\frac{D}{12}}`$, where $`D`$ is the dimensionality of the system). Constant volume is chosen so that $`\rho _{poly}=\frac{\frac{4}{3}_{i=1}^N\pi r^3}{L^2}=1`$. For all of the systems, up to $`1\times 10^8`$ steps of equilibration and $`1\times 10^9`$ steps of data collection are used, providing well converged data. Data accuracy and convergence properties were established through numerous shorter runs.

In several cases presented below it is useful to have an estimate of the location of the putative Kauzmann temperature as calculated through the random first order theory. Here we follow Parisi and coworkers ParisiThermoParisiThermoLJ ; ParisiThermoSoft . Harmonic solid entropies are calculated for $`N`$ particles in $`D`$ dimensions as

$$S_{sol}^{(a)}=\frac{ND}{2}\left(1+log\left(\frac{2\pi }{\beta }\right)\frac{1}{N_{pos}}\underset{i=1}{\overset{N_{pos}}{}}log(|\lambda _i|)\right),$$
(3)
where $`N_{pos}`$ denotes the number of positive eigenvalues $`\lambda _i`$ of the instantaneous Hessian. These are averaged over 100-200 different configurations for each temperature, each separated by $`1\times 10^4`$ steps, which is more than enough to obtain good statistics. Liquid entropies are calculated by thermodynamic integration as

$$S_{liq}=S_{liq}^0+\beta E_{liq}(\beta )_0^\beta 𝑑\beta ^{}E_{liq}(\beta ^{}),$$
(4)
where $`S_{liq}^0`$ represents the ideal gas entropy in the $`\beta 0`$ limit and $`E_{liq}(\beta )`$ is the average energy at inverse temperature $`\beta `$. Liquid entropies are extrapolated into the glassy regime by fitting the data as $`S_{liq}(T)=aT^{2/5}+b`$, a well known fitting form for the entropy of simple liquids Yakov ; ParisiThermoSoft . Inherent structures, when needed, are found by steepest descent quenches for the binary soft-sphere system, and a combination of steepest descent and conjugate gradient for the polydisperse system. Data is averaged over thousands of such configurations for each temperature.

## III Binary Soft-Sphere System

The first system that we investigate is the binary soft-sphere liquid defined in Sec. 2. To demonstrate the robustness of estimates of the Kauzmann temperature in this system, we calculate $`T_K`$ via the random first order theory of Mezard and Parisi ParisiThermoParisiThermoLJ ; ParisiThermoSoft  and the inherent structure-based method of Buchner and Heuer Buchner1 ; Gaussian . The approach of Buchner and Heuer is similar to that of Sciortino, Kob and Tartaglia Kob1 ; Kob2 , and yields nearly the same estimate of $`T_K`$. Using the approach of Mezard and Parisi and the definitions of the entropy of the liquid and disordered solid discussed in Sec. 2, we find $`\mathrm{\Gamma }_K=1.72`$ ($`\mathrm{\Gamma }=\rho T^{\frac{D}{12}}`$) while the method of Buchner and Heuer yields $`\mathrm{\Gamma }_K=1.69`$. A previous estimate using the random first order approach by Coluzzi et. al. ParisiThermoSoft  yielded $`\mathrm{\Gamma }_K=1.65`$ analytically and $`\mathrm{\Gamma }_K=1.75`$ from simulation data. It is unclear why our estimate differs from that of Coluzzi et.al., but it is most likely due to slight differences in numerical implementation. Regardless of these differences, a consistent estimate of $`\mathrm{\Gamma }_K1.7`$ emerges from these approaches.

Recently, Grigera and Parisi implemented the SMC procedure for the soft-sphere mixture. While for large system sizes they could not equilibrate the system at low temperatures (large $`\mathrm{\Gamma }`$), they found a peak in the specific heat of a small (34 particle) system at $`\mathrm{\Gamma }1.7`$. We have reproduced this calculation in Fig. 1a. The agreement between the values of the specific heat calculated from fluctuations and derivatives of the average energy suggest that the equilibration of configurations has been achieved. In the inset of Fig. 1, we display an expanded view of the specific heat for the 34 particle system that includes higher temperatures (lower $`\mathrm{\Gamma }`$ values). A second interesting signature in the specific heat appears near $`\mathrm{\Gamma }=1.25`$. This signature becomes sharper for larger systems, as shown in Fig. 1B. The sharpness of this peak, together with its system-size dependence, strongly support the view that this feature is thermodynamic in nature. It is interesting to note that $`\mathrm{\Gamma }=1.25`$ is close to the onset of supercooling at $`\mathrm{\Gamma }_01.3`$ BrumerPrePrint . For systems that are larger still ($`N>70`$) the signature near $`\mathrm{\Gamma }=1.25`$ vanishes, as does the ability to equilibrate the system with the SMC method for $`\mathrm{\Gamma }>1.5`$ for reasonable trajectory lengths.

The appearance of secondary features in thermodynamic quantities suggests that perhaps the specific heat peak near $`\mathrm{\Gamma }=1.7`$ is not the result of the type of entropy crisis envisioned in the random first order theory of glassy thermodynamics. To investigate this further, we examine the structures obtained from SMC in small systems that display a specific heat peak near $`\mathrm{\Gamma }=1.7`$. In Fig. 2 we show representative configurations of the system obtained by SMC at $`\mathrm{\Gamma }=1.3`$ and $`1.8`$. The system is clearly not an amorphous solid, but a phase separated crystal with substitutional defects, and, at very low temperatures, the phase separated crystal itself.

The ground state structure found by SMC is a phase-separated crystal. By examining the inherent structure configurations sampled at various values of $`\mathrm{\Gamma }`$, we conclude that the peak in the specific heat results from a gap in the density of states between the set of all defective crystal configurations and the ground state of the system. As the system size increases, defective crystalline configurations become harder to locate relative to the set of amorphous configurations, and the sharp features in the specific heat vanish, as does the ability of SMC to equilibrate the system significantly below the estimated value of $`\mathrm{\Gamma }_c`$, the location of the mode-coupling temperature (density)Swapping2 . In Fig. 3 we plot the specific heat for 70 particles obtained from the SMC method. Note that the agreement between the heat capacity calculated via fluctuations and by direct differentiation of the average energy diverge from each other near $`\mathrm{\Gamma }_c1.45`$.

We thus conclude that the specific heat peak found in the binary soft-sphere system is not a direct consequence of an ideal glass transition. One may speculate that the gap between defective crystalline states and the ground state occurs near $`\mathrm{\Gamma }_K=1.7`$ because the number of amorphous configurations at these energy values (nearly) vanish, and thus the peak at $`\mathrm{\Gamma }_K=1.7`$ is indeed an indirect consequence of the vanishing of the configurational entropy associated with glassy states. The fact that slightly larger systems (i.e. $`N>70`$) than those for which converged specific heats may be obtained cannot be equilibrated significantly below $`\mathrm{\Gamma }_c\mathrm{\Gamma }_K`$ makes this reasoning suspect. It still remains to be determined if the fact that the predicted Kauzmann value $`\mathrm{\Gamma }_c1.7`$ coincides with the gap between defective crystal states and the ground state is purely accidental. Furthermore, we find that, while SMC is more efficient than methods like parallel tempering for $`\mathrm{\Gamma }\mathrm{\Gamma }_c`$, neither approach is able to equilibrate large systems significantly below the mode-coupling temperature (density) in this system.

## IV Polydisperse Systems

Having seen that states with crystalline order interfere with a direct investigation of the entropy crisis in the binary soft-sphere system, we now turn to the analysis of polydisperse systems. The first goal will be to re-investigate the completely polydisperse two dimensional hard-sphere system studied by Santen and KrauthKrauth . Here we will examine systems in three dimensions, as well as systems at finite temperatures to address concerns that might arise with Santen and Krauth’s particular system.

First, we calculate the compressibility of a 256 particle two dimensional hard-sphere system via the SSMC approach. In Fig. 4 we show $`\kappa ^1`$ vs. $`\rho `$ for this system. Our results agree with those of Santen and Krauth, although we were able to equilibrate the system at slightly higher densities. Interestingly, we have found that our implementation of SSMC is more efficient than the cluster algorithm of Santen and Krauth, although this may indeed be a result of our own inefficient optimization of their technique. Regardless, the similarity in the convergence of the two methods is significant, as the SSMC approach works in arbitrary dimensions, while the method of Santen and Krauth is difficult to generalize, particularly to systems with $`D>2`$ KrauthPersonal . By fitting the diffusion constant of the larger discs to a power law form, $`D(\rho )(\rho \rho _c)^\gamma ,`$ Santen and Krauth extract a “glass transition” density of $`\rho _c0.805`$. These diffusion constants are measured by short-time local Monte Carlo dynamics, which has been shown to accurately yield relaxation times in two dimensional systems Harrowell  and meaningful quantities in the $`\alpha `$-relaxation regime Gleim . This means of extracting a critical density should actually yield the mode-coupling density and not the glass transition density. Thus, it is not surprising that there is no thermodynamic signature in the compressibility near $`\rho _c=0.805`$. One can make a rough estimate of the location of $`\rho _K`$ in the system of Santen and Krauth by noting that a conservative guess for the ratio of mode-coupling to Kauzmann densities from the $`\rho T^{\frac{D}{12}}`$ dependence in an “equivalent” soft-sphere system would imply that $`\rho _K`$ is about 10 percent larger than $`\rho _c`$. While this estimate is crude, it implies that exploration of the putative entropy crisis in the system studied by Santen and Krauth is out of reach of both the SSMC and cluster algorithm used by Santen and Krauth.

To further explore the properties of this two dimensional hard-sphere system, we address the nature of the phase space available to the smaller discs. While, at the highest densities we can probe, the largest discs are completely trapped for long Monte Carlo times, the smaller discs are diffusive, and do not even show a significant cage effect. In Fig. 5 we show the mean square displacement of large and small discs, illustrating that part of the system remains fluid even at the highest densities. From the standpoint of the standard definition of configurational entropy, it would appear that the polydisperse hard-sphere system studied by Santen and Krauth cannot experience a strict entropy crisis in principle due to the disparity of disc sizes. Recent work by Kumar et. al. suggests that a reasonable correlation should exist between the density of random close-packing and the Vogel density in a hard sphere system Kumar . If we crudely interchange the Vogel density and the Kauzmann density, then we expect that the ideal glass transition, if one occurs at all, will occur close to the random close-packed density. If, as in the polydisperse system of Santen and Krauth, the larger spheres freeze while the smaller spheres remain mobile, then the existence of the mobile, incompressible background fluid will always prevent an effective close-packing of the larger spheres. It is not clear what signatures of an entropy crisis might exist in systems where a finite fraction of particles may still undergo diffusive motion while the remaining particles are frozen, but it would appear that the very concept of an entropy crisis in such a system is perhaps ill-defined.

Since the system studied by Santen and Krauth is both athermal and two dimensional, we have also studied a three dimensional soft-sphere version of this system. In Fig. 6 we show the specific heat of a small (N=34) version of the polydisperse three dimensional soft-sphere system discussed in the introduction. Interestingly, while the data shows that the heat capacity calculated via fluctuations and by derivatives of the average energy agree, there is a sharp thermodynamic feature in the heat capacity at low temperatures. One would naively expect that this system cannot exhibit a state with intermediate range order and therefore this peak must be related to a glassy phenomenon, but this is not the case. Subtle features in the radial distribution function as well as close inspection of instantaneous and inherent structures make it clear that particles of similar size have bunched together. In Fig. 7 we show one such inherent configuration found via SSMC at low temperatures.

We conclude once again that the peak in the heat capacity is the result of an exotic crystallization phenomena. In this sense, our three dimensional thermalized version of the system studied by Santen and Krauth is similar to the binary soft-sphere system studied by Grigera and Parisi. For larger systems the range of thermodynamic parameters for which non-local Monte Carlo methods provide equilibrated results is not broad enough to investigate deeply glassy states, while smaller systems show thermodynamic signatures consistent with the formation of defective or phase separated crystals. Interestingly, we could not find evidence of this type of phase separation in small, two dimensional polydisperse hard or soft-sphere systems, however we were still limited to densities and temperatures that differ significantly from those at which an entropy crisis might be expected.

## V Conclusions

In this work we have taken a closer look at the recent use of non-local Monte Carlo methods to investigate the nature of low temperature glassy thermodynamics. We have found that claims for and against the existence of an entropy crisis that have been recently reported in the literature must be reconsidered. The peak in the specific heat found by Grigera and Parisi through the SMC approach can be attributed to an energy gap between the manifold of defective crystalline configurations and the phase separated ground state of the system. The use of SSMC (optimized SMC) performs at least as well as the cluster algorithm of Santen and Krauth for the two dimensional polydisperse hard-sphere system, but both are unable to probe this system at high enough densities to render strong claims about the behavior of the configurational entropy meaningful. Furthermore, since this system has a very wide range of particle radii, a macroscopic portion of the system remains fluid at the highest densities. This fact further clouds the discussion of a possible entropy crisis in this system. It is not surprising that even powerful Monte Carlo methods that make use of non-local moves still cannot equilibrate systems under conditions where the number of available configurations become sparse. Interestingly, while the Monte Carlo methods investigated here greatly increase the rate of equilibration for temperatures (densities) near $`T_c`$ ($`\rho _c`$) they are not efficient at lower temperatures (higher densities). Completely polydisperse systems are an exception to this finding, but these systems are pathological in the sense that there is always significant phase space for a finite fraction of the system to explore even when part of the system is configurationally frozen. In general, we have found that one has to be careful to avoid the subtle phase separation and the corresponding sharp features in thermodynamic quantities that may occur in these systems.

While the findings presented here are somewhat negative, they do provide important guidelines and warning signs for any computational investigation of the entropy crisis in systems with realistic potentials. It would be interesting to perform non-local Monte Carlo studies on small systems whose potentials are tailored to strongly penalize structures with intermediate and long range order Leonardo . It would also be interesting to re-investigate larger versions of the systems studied here with the more powerful Wang-Landau approach WangWang2 ; Faller . We are currently pursuing these directions.

###### Acknowledgements.

We gratefully acknowledge NSF support through grant #CHE-0134969 and an NSF fellowship to Y.B.. D.R.R. is and Alfred P. Sloan Foundation Fellow and a Camille Dreyfus Teacher-Scholar.

Figure Captions

Fig. 1: (a) The heat capacity per particle for a 34 particle binary soft-sphere system. The inset extends this calculation to higher temperatures (lower $`\mathrm{\Gamma }`$), where a second smaller peak is evident. $`C_V`$ as calculated by fluctuations is plotted as (x) and finite difference as (\*), a convention which is maintained throughout the paper. Swap moves were attempted on average every ten moves, and data was collected for between $`1\times 10^7`$and $`1\times 10^9`$ steps, with longer runs for lower temperatures. (b) The heat capacity per particle for the same system and computational details as in part (a), but with 60 particles. The heat capacity as calculated by fluctuations and derivatives is shown. The 34 particle peak is included (solid line) to facilitate direct comparison.

Fig. 2: Crystalline structures found by SMC. (a) shows an inherent configuration of a 60 particle binary soft- sphere system sampled at $`\mathrm{\Gamma }=1.3`$, near the first heat capacity peak and (b) presents an instantaneous structure sampled at $`\mathrm{\Gamma }=1.8`$, near the second peak.

Fig. 3: The heat capacity per particle for the 70 particle binary soft-sphere system. Equilibrated samples are obtained only up to $`\mathrm{\Gamma }=1.45`$. Swap moves were attempted on average every ten moves, and data was collected for between $`1\times 10^8`$ and $`1\times 10^9`$ steps.

Fig. 4: The converged inverse isothermal compressibility $`\kappa ^1`$ vs. $`\rho `$ for the polydisperse hard-sphere system of Santen and Krauth. The system was equilibrated by SSMC. Swap moves were attempted on average every ten moves and data was collected for $`1\times 10^9`$ steps, though such a long run was only necessary for extremely high pressures.

Fig. 5: Mean square displacement vs. MC time for two small and two large particles chosen at random from the polydisperse hard-sphere system of Santen and Krauth. The data is taken at $`P=0.09`$, which corresponds to $`\rho 0.861`$.

Fig. 6: The heat capacity per particle for a 34 particle polydisperse soft-sphere system. The peak arises due to a phase separation process that induces intermediate range order in the system. Swap moves are attempted on average every ten steps and data was collected for $`1\times 10^7`$ to $`6\times 10^8`$ steps.

Fig. 7: An inherent structure of the 34 particle polydisperse soft-sphere system at $`\mathrm{\Gamma }=10.0`$. A cluster of large particles lies in the foreground, and a cluster of small particles inside (not behind) the large cluster. Similar clustering and the beginnings of phase separation can be found in all low temperature structures of this system.

Fig. 1a

Fig. 1b

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7