format: text

text:
# Measuring velocity of sound with nuclear resonant inelastic x-ray scattering

## Appendix

Our study of long-wavelength vibrational excitations starts with the continuum version of the equations of motion of a set of harmonically bound atoms, which reads

$$\omega ^2m^\mu 𝐮^\mu =\underset{\nu }{}\mathrm{\Phi }^{\mu \nu }𝐮^\nu ,$$
(15)
where $`\mathrm{\Phi }^{\mu \nu }`$ are the force constant matrices and $`𝐮^\mu `$ is the displacement vector of atom $`\mu `$ for a vibrational mode of energy $`\mathrm{}\omega `$. We introduce mass and atomic density functions $`\rho (𝐱)`$ and $`\eta (𝐱)`$ by

$`\rho (𝐱)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}m^\mu \delta ^3(𝐱𝐫_\mu )`$
$`\eta (𝐱)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\delta ^3(𝐱𝐫_\mu ),`$ (16)
where $`𝐫_\mu `$ is the position of atom $`\mu `$. The continuum version that substitutes Eq. (15) is then given by

$$\omega ^2\rho (𝐱)𝐮(𝐱)=\eta (𝐱)\mathrm{\Phi }(𝐱𝐱^{})\eta (𝐱^{})𝐮(𝐱^{})d^3x^{},$$
(17)
where we made the usual assumption that the force constant matrix depends on coordinate differences only. The original intentions to study long-wavelength excitations are served best by introducing Fourier transforms, e.g., $`\stackrel{~}{𝐮}(𝐪)=FT[𝐮(𝐱)]`$, and eventually expanding for small values of momentum. The transformed Eq. (17) reads

$`\omega ^2`$ $`{\displaystyle }\stackrel{~}{\rho }(𝐪𝐪^{})\stackrel{~}{𝐮}(𝐪^{})d^3q^{}={\displaystyle \frac{1}{(2\pi )^3}}\times `$ (18)
$`{\displaystyle \stackrel{~}{\eta }(𝐪𝐪^{})\stackrel{~}{\mathrm{\Phi }}(𝐪^{})\stackrel{~}{\eta }(𝐪^{}𝐪^{\prime \prime })\stackrel{~}{𝐮}(𝐪^{\prime \prime })d^3q^{}d^3q^{\prime \prime }}.`$
The quantities $`\stackrel{~}{\eta }(𝐪)`$ and $`\stackrel{~}{\rho }(𝐪)`$ are closely related to the structure function $`S(𝐪)`$ that is typically obtained from x-ray or neutron diffraction experiments. For crystals, these functions are described by a series of very sharp peaks at values given by the reciprocal lattice vectors. Disordered or amorphous materials do not produce these sharp peaks with the exception of the $`𝐪=0`$ maximum which is not related to spatial order. One can assume that some of the salient features of Eq. (18) are captured by retaining the $`𝐪=0`$ maximum only, i.e., we approximate

$`\stackrel{~}{\rho }(𝐪)`$ $`=`$ $`(2\pi )^3{\displaystyle \frac{m}{V}}\delta ^3(𝐪)`$
$`\stackrel{~}{\eta }(𝐪)`$ $`=`$ $`(2\pi )^3{\displaystyle \frac{1}{V}}\delta ^3(𝐪),`$ (19)
where $`m`$ and $`V`$ are the average mass and volume per atom. Equation (18) simplifies to

$$m\omega ^2\stackrel{~}{𝐮}(𝐪)=\frac{1}{V}\stackrel{~}{\mathrm{\Phi }}(𝐪)\stackrel{~}{𝐮}(𝐪),$$
(20)
and describes a class of solutions that should be common for all harmonic solids. The previous equation permits us to label the modes with the value of $`𝐪`$ and a branch index $`s`$ originating from the tensor character of the force-constant matrix, i.e., $`\omega =\omega _{𝐪s}`$. We note that $`\stackrel{~}{\mathrm{\Phi }}(0)=0`$ follows from the invariance of the vibrational energy with respect to a displacement identical for all atoms. Assuming that the force-constant matrix falls off sufficiently fast with distance an expansion of Eq. (20) in powers of $`q=|𝐪|`$ is permissible. The first-order term in this expansion will also vanish under the reasonable assumption of inversion symmetry of the force-constant matrix. For hydrodynamic modes, i.e., $`q0`$, we may therefore use

$$\stackrel{~}{\mathrm{\Phi }}(𝐪)=\frac{q^2}{2}\left(\frac{^2\stackrel{~}{\mathrm{\Phi }}}{q^2}\right)_{q=0},$$
(21)
and long-wavelength vibrational excitations are described by a modified Eq. (20)

$$m\omega _{𝐪s}^2\stackrel{~}{𝐮}(𝐪)=q^2\stackrel{~}{\mathrm{\Theta }}(\widehat{𝐪})\stackrel{~}{𝐮}(𝐪).$$
(22)
We obtain linearly dispersing modes with energies $`\mathrm{}\omega _{𝐪s}=c_{\widehat{𝐪}s}\mathrm{}q`$ with sound velocities $`c_{\widehat{𝐪}s}`$ that depend on the direction of $`𝐪`$. The matrix $`\stackrel{~}{\mathrm{\Theta }}`$ is related to the elastic tensor of the material (see for example ref. Ashcroft and Mermin, 1976) and is given by

$$\stackrel{~}{\mathrm{\Theta }}(\widehat{𝐪})=\frac{1}{2V}\left(\frac{^2\stackrel{~}{\mathrm{\Phi }}}{q^2}\right)_{q=0}=\frac{1}{2N}\underset{\mu }{}(\widehat{𝐪}𝐫_\mu )^2\mathrm{\Phi }(𝐫_\mu ).$$
(23)
$`N`$ is the number of atoms in the sample. The eigenvalues of $`\stackrel{~}{\mathrm{\Theta }}`$ provide us with the sound velocities, and the corresponding eigenvectors will not depend on $`q`$. The displacement field for a mode $`𝐪s`$ is then given by

$$𝐮_{𝐪s}^\mu =\alpha _{𝐪s}𝐩_{\widehat{𝐪}s}e^{i𝐪𝐫_\mu },$$
(24)
where $`\widehat{𝐩}_{\widehat{𝐪}s}`$ is a normalized eigenvector of $`\stackrel{~}{\mathrm{\Theta }}`$ describing the polarization of the mode, and $`\alpha _{𝐪s}`$ is an appropriately chosen normalization factor. Equation (24) shows that, even for disordered materials, hydrodynamic modes are equivalent to plane-wave excitations of the atomic displacements.


text:
# Lower Order Terms in Szegö Theorems on Zoll Manifolds

## 1. Introduction and main results

The main motivation for this work was to find an explicit formula for a “Szegö–regularized” determinant of a zeroth order pseudodifferential operator (PsDO) on a Zoll manifold introduced in \[GO1, O2\], see Remark 2.4. Our main result, Theorem 1.2, is valid for any dimension $`d`$. In the case $`d=2`$, Theorem 1.2 gives such a formula.

In this paper we find a third order generalization of the Strong Szegö Limit Theorem (SSLT) for a zeroth order PsDO on the unit circle (Theorem 1.1), and on a Zoll manifold of an arbitrary dimension $`d`$ (Theorem 1.2). We give also an outline of the proof of Theorem 1.2 and derive Theorem 1.1 from the former. The detailed proof of Theorem 1.2 can be found in \[Gi1, Chapter 1\] and will be published in a forthcoming paper \[Gi2\]. In Section 3, we sketch a proof of the generalized Hunt–Dyson combinatorial formula (Theorem 3.1), see \[Gi3\] for the detailed proof. We give also a brief review of the related work.

Recall that $`M`$ is called a Zoll manifold if it is compact, closed and such that the geodesic flow on $`M`$ is simply periodic with period $`2\pi `$. The unit circle and the standard sphere of any dimension are Zoll manifolds. A second order generalization of the SSLT for a Zoll manifold of any dimension has been obtained by V. Guillemin and K. Okikiolu \[GO1, GO2\], see also an important preceding work \[O1\] by K. Okikiolu for the case of the two- and three-dimensional sphere. The proofs in \[O1, GO1, GO2\] use a combinatorial identity due to G. A. Hunt and F. J. Dyson and proceed in the spirit of the combinatorial proof of the classical SSLT by M. Kac \[K\]. See also \[GO3, O2\] where the combinatorial approach and the Hunt–Dyson formula are used in a different setting to obtain a second order generalization of the SSLT for a manifold with the set of closed geodesics of measure zero in the unit cotangent bundle.

In the proof of Theorem 1.2 we use the method of \[GO2\]. A central role in our proof is played by a certain combinatorial identity which generalizes the Hunt–Dyson formula mentioned above to an arbitrary natural power. We call this identity the generalized Hunt–Dyson formula (gHD), see Theorem 3.1.

After having discovered and proved the formula gHD we realized that it is related to another combinatorial theorem, which has a long history. The mentioned theorem is a result due to H. F. Bohnenblust that appeared in an article by F. Spitzer on random walks \[S, Theorem 2.2\], and is now commonly known as the Bohnenblust–Spitzer theorem (BSt). The characteristic function of the maximum of a random walk with independent identically distributed steps has been computed in \[S\] with the help of the BSt. On the other hand, using the gHD, we can compute \[Gi3\] the moment of an arbitrary order of the mentioned maximum (note that the usual Hunt–Dyson formula allows one to compute only the expected value of the maximum, see \[K\]). This indicates that the gHD and the BSt should be essentially the same. And indeed this is the case: it turns out that the gHD can be derived from the BSt, and vice versa, see \[Gi3\].

We have found the formula gHD being unaware of the BSt. In \[Gi1, Chapter 2\], a proof of the gHD “from scratch” can be found. This together with a a derivation of the BSt from the the gHD in \[Gi1, Gi3\] gives a new proof of the BSt. There are known at least four other proofs: the original one \[S, F\], a proof by G. Baxter \[B1, B2\], by J. G. Wendel \[We\], and finally, a unifying approach of G.-C. Rota in the framework of (Glen) Baxter algebras \[Ro, RoSm\]. See Section 3 for more combinatorial details.

Let $`det`$ denote the determinant of a finite rank operator. In Section 2, we give explicit asymptotic formulas for $`\mathrm{log}detP_nBP_n`$, $`n\mathrm{}`$, where $`P_n`$ is the projection onto the first $`n`$ eigenspaces of the Laplace–Beltrami operator on a Zoll manifold, and $`B`$ is an arbitrary zeroth order PsDO, see Corollary 2.2 (for $`d=1`$) and Corollary 2.3 ($`d2`$).

Let $`𝕊^1`$ be the unit circle $`/2\pi `$. Denote by $`P_n`$, $`n`$, the orthogonal projection from $`L^2(𝕊^1)`$ to the subspace spanned by $`\{e^{ikx}\}_{|k|n}`$. For a function $`fL^1(𝕊^1)`$ denote its $`k`$th Fourier coefficient by

(1.1)
$$\widehat{f}_k:=_0^{2\pi }f(x)e^{ikx}\frac{dx}{2\pi },k.$$

Let $`b(x)`$ be a strictly positive function on $`𝕊^1`$ such that $`_k|k||\widehat{(\mathrm{log}b)}_k|^2<\mathrm{}`$. Denote by $`B`$ the operator of multiplication by $`b`$ acting in $`L^2(𝕊^1)`$. Recall that the matrix representation of the operator $`B`$ is the Toeplitz matrix $`(\widehat{b}_{jk})_{j,k}`$. The classical Strong Szegö Limit Theorem (SSLT) \[Sz\] states that

$$Tr\mathrm{log}P_nBP_n=TrP_n(\mathrm{log}B)P_n+\underset{k=1}{\overset{\mathrm{}}{}}k\widehat{(\mathrm{log}b)}_k\widehat{(\mathrm{log}b)}_k+o(1),n\mathrm{}.$$

Here $`Tr\mathrm{log}P_nBP_n=\mathrm{log}detP_nBP_n`$ and $`TrP_n(\mathrm{log}B)P_n=(2n+1)_0^{2\pi }\mathrm{log}b(x)\frac{dx}{2\pi }`$. It has been shown by H. Widom that the remainder is $`O(n^{\mathrm{}})`$ if $`b(x)C^{\mathrm{}}(𝕊^1)`$, see \[W1\].

Let $`M`$ be a closed manifold of dimension $`d`$. Let $`\mathrm{\Psi }^m(M)`$, $`m`$, denote the space of the classical PsDO of order $`m`$ on $`M`$. Recall that for a given $`G\mathrm{\Psi }^m(M)`$, its principal symbol $`\sigma _m(G)`$ and subprincipal symbol $`sub(G)`$ are well-defined on $`T^{}M`$. Let $`S^{}M`$ be the unit cotangent bundle and denote by $`dxd\xi `$ the standard measure on $`S^{}M`$ divided by $`(2\pi )^d`$.

The simplest form of our result is for $`d=1`$ and the function $`\psi (u)=\mathrm{log}u`$ which is analytic in a neighborhood of $`u=1`$.

###### Theorem 1.1.

Let $`M=𝕊^1`$ and $`P_n`$ be the projection on the linear span of $`\{e^{ikx}\}_{|k|n}`$. Let $`B\mathrm{\Psi }^0(M)`$ and assume that $`\sigma _0(B)`$ is strictly positive and a certain symbolic norm of $`IB`$ is sufficiently small. Then $`\mathrm{log}B\mathrm{\Psi }^0(M)`$ and the following holds

(1.2) $`Tr\mathrm{log}`$ $`P_nBP_n=TrP_n(\mathrm{log}B)P_n+{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\sigma _0(\mathrm{log}B))}_k\widehat{(\sigma _0(\mathrm{log}B))}_kdxd\xi `$
$`+{\displaystyle \frac{1}{n}}{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\sigma _0(\mathrm{log}B))}_k\widehat{(sub(\mathrm{log}B))}_kdxd\xi +O\left({\displaystyle \frac{1}{n^2}}\right),n\mathrm{}.`$
In (1.2) the argument $`(x,\xi )S^{}M`$ is omitted for brevity. By the Fourier coefficient we mean the following: for a fixed $`(x,\xi )S^{}M`$ compute the Fourier coefficient in accordance with (1.1) along the unit circle being the closed geodesic starting at $`(x,\xi )`$.

Let us denote $`b_0:=\sigma _0(B)`$ and $`b_{\mathrm{sub}}:=sub(B)`$. Recall that for an analytic $`f`$, $`\sigma _0(f(B))=f(\sigma _0(B))`$ and $`subf(B)=f^{}(\sigma _0(B))sub(B)`$. Then (1.2) takes the form

(1.3) $`Tr\mathrm{log}P_nBP_n`$ $`=TrP_n(\mathrm{log}B)P_n+{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\widehat{(\mathrm{log}b_0)}_kdxd\xi `$
$`+{\displaystyle \frac{1}{n}}{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\widehat{(b_{\mathrm{sub}}/b_0)}_kdxd\xi +O\left({\displaystyle \frac{1}{n^2}}\right),n\mathrm{}.`$
###### Remark 1.1.

Theorem 1.1 is a particular case of Theorem 1.2 below, see Remark 1.4.

Let us fix some notations and then state the result for the higher dimensional case. Let $`M`$ be a Zoll manifold of dimension $`d`$. Let $`\mathrm{\Delta }`$ denote the Laplace–Beltrami operator on $`M`$. It is known \[DG\] that there exists a constant $`\alpha `$ such that the spectrum of $`\sqrt{\mathrm{\Delta }}`$ lies in bands around the points $`k+\frac{\alpha }{4}`$, $`k`$. Moreover, in has been shown in \[CdV\] that there exists $`A_1\mathrm{\Psi }^1(M)`$ such that $`[\mathrm{\Delta },A_1]=0`$ and the spectrum of the operator

(1.4)
$$A:=\sqrt{\mathrm{\Delta }}\frac{\alpha }{4}A_1$$

is $``$. Let $`P_n`$, $`n`$, denote the projection from $`L^2(M)`$ onto the subspace spanned by the eigenfunctions of $`A`$ corresponding to the eigenvalues $`1,2,\mathrm{},n`$. Let $`dxd\xi `$ be the defined above measure on $`S^{}M:=\{(x,\xi ):\sigma _1(A)(x,\xi )=1\}`$. Following \[GO1, GO2\], for a $`B\mathrm{\Psi }^0(M)`$ and $`t`$ introduce the operator

(1.5)
$$B^t:=e^{itA}Be^{itA}.$$

By Egorov’s theorem, $`B^t\mathrm{\Psi }^0(M)`$, and also

$$\sigma _0(B^t)(x,\xi )=\sigma _0(B)(\mathrm{\Theta }^t(x,\xi )),$$

where $`\mathrm{\Theta }^t`$ stands for the shift by $`t`$ units along geodesic flow. Note that because $`sub(A)=const`$ the following also holds \[G, Lemma 2.2\]

$$sub(B^t)(x,\xi )=sub(B)(\mathrm{\Theta }^t(x,\xi )).$$

Because $`spec(A)=`$, the operator $`B^t`$ is periodic in $`t`$ with period $`2\pi `$.

The result we are about to state, and also Theorem 1.1 above, gives information on the asymptotic behavior of

(1.6)
$$Tr\psi (P_nBP_n)TrP_n\psi (B)P_n,n\mathrm{},$$

for the analytic in a neighborhood of $`1`$ function $`\psi (u)=\mathrm{log}u`$. In \[Gi1, Gi2\] we obtain the corresponding statements for an arbitrary analytic function $`\psi (u)`$. These results involve certain maps $`W`$, $`\mathrm{\Phi }`$, $`\mathrm{{\rm Y}}`$ whose action on an analytic function is a continuous function from $`^j`$ to $``$ for some $`j`$. The map

$$W[\psi ](y_1,y_2):=\frac{1}{2}_{y_1}^1_{y_2}^1\frac{\psi ^{}(u_1)\psi ^{}(u_2)}{u_1u_2}𝑑u_1𝑑u_2,y_1,y_2>0,$$

was obtained by A. Laptev, D. Robert and Yu. Safarov in \[LRS\] (here $`j=2`$). Earlier, an equivalent to $`W`$ map has been obtained by H. Widom \[W2, W3\] in a two-term Szegö type asymptotics for integral operators with discontinuous symbols. It has been noticed by H. Widom, and also by the authors \[LRS, O2\], that for the function $`\psi (u)=\mathrm{log}u`$, the action $`W[\mathrm{log}]`$ takes a simple form

$$W[\mathrm{log}](y_1,y_2)=\frac{1}{2}\mathrm{log}y_1\mathrm{log}y_2.$$

This explains the fact that the map $`W`$ is not required in the next theorem.

Now define for $`y_1,y_2,y_3<1`$,

(1.7) $`\mathrm{\Phi }[\mathrm{log}](y_1,y_2,y_3)`$ $`=y_3{\displaystyle _0^{y_1}}{\displaystyle _0^{y_2}}[{\displaystyle \frac{\mathrm{log}(1u_1)}{u_1(u_1u_2)(u_1u_3)}}`$
$`{\displaystyle \frac{\mathrm{log}(1u_2)}{u_2(u_1u_2)(u_2u_3)}}+{\displaystyle \frac{\mathrm{log}(1y_3)}{y_3(u_1y_3)(u_2y_3)}}]du_1du_2.`$
For the map $`\mathrm{\Phi }`$, $`j=3`$. The expression for $`\mathrm{{\rm Y}}[\mathrm{log}]`$ appearing in the statement of the next theorem is complicated and is not given here, see \[Gi1, (1.1.31)\] and \[Gi2\] for details. It depends only on the principal symbol $`\sigma _0(B)`$ and involves certain Poisson brackets of the type $`\{\sigma _0(B^t),\sigma _0(B^s)\}`$, $`0t,s2\pi `$.

For a function $`fC^{\mathrm{}}(S^{}M)`$ introduce the $`k`$th Fourier coefficient along the closed geodesic of length $`2\pi `$ starting at a given point $`(x,\xi )`$

$$\widehat{f}_k(x,\xi ):=_0^{2\pi }e^{ikt}f(\mathrm{\Theta }^t(x,\xi ))\frac{dt}{2\pi },k.$$

For an arbitrary $`B\mathrm{\Psi }^0(M)`$, let us write $`b_0:=\sigma _0(B)`$, $`b_{\mathrm{sub}}:=sub(B)`$, $`b_0^t(x,\xi ):=b_0(\mathrm{\Theta }^t(x,\xi ))`$, and omit the argument $`(x,\xi )S^{}M`$ for brevity.

###### Theorem 1.2.

Let $`M`$ be a Zoll manifold of dimension $`d`$. Let $`A`$ be defined by (1.4). Assume that $`\sigma _1(A)(x,\xi )=\sigma _1(A)(x,\xi )`$ for all $`(x,\xi )T^{}M`$. Let $`B\mathrm{\Psi }^0(M)`$ be such that $`\sigma _0(B)`$ is strictly positive and a certain symbolic norm of $`IB`$ is sufficiently small. Then $`\mathrm{log}B\mathrm{\Psi }^0(M)`$ and the following holds

(1.8) $`Tr\mathrm{log}`$ $`P_nBP_n=TrP_n(\mathrm{log}B)P_n`$
$`+n^{d1}{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\widehat{(\mathrm{log}b_0)}_kdxd\xi `$
$`+n^{d2}{\displaystyle \frac{1}{2}}{\displaystyle _{S^{}M}}({\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\widehat{(b_{\mathrm{sub}}/b_0)}_k`$
$`+(d1)[{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}(k^2+(1+\alpha /2)k)\widehat{(\mathrm{log}b_0)}_k\widehat{(\mathrm{log}b_0)}_k`$
$`+{\displaystyle \underset{k,l=1}{\overset{\mathrm{}}{}}}kl{\displaystyle _0^{2\pi }}{\displaystyle _0^{2\pi }}{\displaystyle _0^{2\pi }}e^{ik(tr)+il(sr)}\mathrm{\Phi }[\mathrm{log}](b_0^t,b_0^s,b_0^r){\displaystyle \frac{dt}{2\pi }}{\displaystyle \frac{ds}{2\pi }}{\displaystyle \frac{dr}{2\pi }}]`$
$`+\mathrm{{\rm Y}}[\mathrm{log}](b_0))dxd\xi `$
$`+O(n^{d3}),n\mathrm{}.`$
###### Remark 1.2.

The existence of a full expansion of the type (1.8) has been proven in \[GO1\]. Explicit formulas for the coefficients of the first first two terms in (1.8) (and also of the $`\mathrm{log}n`$ term which is a part of $`TrP_n(\mathrm{log}B)P_n`$) have been obtained in \[GO1, GO2\].

###### Remark 1.3.

The formula (1.7) for $`\mathrm{\Phi }[\mathrm{log}]`$ has a structure similar to the coefficient in the third asymptotic term in a Szegö type expansion for convolution operators obtained by R. Roccaforte in \[R\]. A version of the Bohnenblust–Spitzer combinatorial theorem is used in the proof in \[R\]. A second order Szegö type expansion for convolution operators was established by H. Widom in \[W5\] with the help of the usual Hunt–Dyson combinatorial formula (3.4). A full asymptotic expansion for convolution operators has been found in \[W6\].

###### Remark 1.4.

For $`d=1`$, the square bracket in (1.8) disappears. Also $`\mathrm{{\rm Y}}`$ vanishes, because all the Poisson brackets vanish in this case (for each of the two cotangent directions the angle does not change and $`\sigma _0(B)`$ is homogeneous of degree $`0`$ in $`\xi `$). Theorem 1.1 follows.

###### Remark 1.5.

As we have noticed above, Theorem 1.1 and Theorem 1.2 give lower order corrections to (1.6) for $`\psi (u)=\mathrm{log}u`$, as $`n\mathrm{}`$. In the next section we find explicit formulas for the coefficients in the asymptotic expansion of $`Tr\mathrm{log}P_nBP_n`$, as $`n\mathrm{}`$. These expansions are fourth order for $`d=1,2,3`$ and third order for $`d4`$, as $`n\mathrm{}`$.

###### Proof of Theorem 1.2: an outline.

We start by expanding $`\psi (u)=\mathrm{log}u`$ in a power series about the point $`1`$. After that it suffices to prove (1.8) for $`\psi (u)=u^m`$ for all $`m`$.

Let us recall the method of \[GO2\]. Let $`\pi _k`$, $`k`$, be the projection on the $`k`$th eigenspace of the operator $`A`$ and set $`\pi _k:=0`$ for $`k0`$. Then $`P_n=_{k=1}^n\pi _k`$ for $`n`$, and we set $`P_n:=0`$, $`n0`$. Recall that because $`M`$ is a Zoll manifold, for an arbitrary $`B\mathrm{\Psi }^0(M)`$, the operator $`B^t`$, $`t`$, defined by (1.5) is $`2\pi `$-periodic. Therefore one can introduce the Fourier expansion $`B=_jB_j`$ where $`B_j\mathrm{\Psi }^0(M)`$, $`j`$, is defined by

$`B_j`$ $`={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\pi _{k+j}B\pi _k`$ $`={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}e^{ijt}e^{itA}Be^{itA}𝑑t.`$
For $`m`$ and $`j_1,\mathrm{},j_m`$ introduce the notation

$$M_m(\overline{j}):=min(0,j_1,j_1+j_2,\mathrm{},j_1+\mathrm{}+j_m).$$

The following commutation relation is of central importance here

(1.9)
$$B_jP_n=P_{n+j}B_j,n,j.$$

Using (1.9) one moves all the projectors to the left in the expression

$$(P_nBP_n)^m=\underset{j_1,\mathrm{},j_m}{}P_nB_{j_1}P_nB_{j_2}P_n\mathrm{}P_nB_{j_m}P_n$$

obtaining $`P_nB^mP_n`$ plus a correction. This implies for all $`n`$

(1.10) $`Tr(P_nB`$ $`P_n)^mTrP_nB^mP_n`$
$`={\displaystyle \underset{j_1+\mathrm{}+j_m=0}{}}Tr\left((P_nP_nP_{n+j_1}\mathrm{}P_{n+j_1+\mathrm{}+j_m})B_{j_1}\mathrm{}B_{j_m}\right)`$
$`={\displaystyle \underset{j_1+\mathrm{}+j_m=0}{}}{\displaystyle \underset{k=0}{\overset{M_m(\overline{j})+1}{}}}Tr(\pi _{n+k}B_{j_1}\mathrm{}B_{j_m}).`$
Next, it is important that for any $`G\mathrm{\Psi }^0(M)`$, $`M`$ being a Zoll manifold, there exists a full asymptotic expansion for $`Tr(\pi _kG)`$, as $`k\mathrm{}`$, see Lemma 1.3 below. This result is due to Y. Colin de Verdière \[CdV\]. The coefficients in that expansion are certain Guillemin–Wodzicki residues. Recall that for any compact closed manifold $`M`$ of dimension $`d`$ the Guillemin–Wodzicki residue of a pseudodifferential operator $`G\mathrm{\Psi }^m(M)`$ of order $`m`$ is defined by

$$Res(G):=_{S^{}M}\sigma _d(G)(x,\xi )𝑑x𝑑\xi .$$

For an arbitrary $`G\mathrm{\Psi }^0(M)`$ denote

$$R_l(G):=Res(A^{d+l}G),l=0,1,2,\mathrm{}.$$

###### Lemma 1.3.

Let $`M`$ be a Zoll manifold of dimension $`d`$. Assume $`G\mathrm{\Psi }^0(M)`$. Then for any $`N=0,1,2,\mathrm{}`$, there exists $`C_N<\mathrm{}`$ such that

(1.11)
$$\left|Tr(\pi _kG)\underset{l=0}{\overset{N}{}}k^{d1l}R_l(G)\right|C_Nk^{d2N},k.$$

See \[CdV\] and \[GO2, Appendix\] for the proof of (1.11).

Now assume that for some $`K<\mathrm{}`$ we have $`B_k=0`$ for $`|k|>K`$ (this is not assumed in the full proof in \[Gi2\]). Then (1.10) and Lemma 1.3 imply

(1.12) $`Tr(`$ $`P_nBP_n)^mTrP_nB^mP_n`$
$`={\displaystyle \underset{j_1+\mathrm{}+j_m=0}{}}\{n^{d1}M_m(\overline{j}){\displaystyle _{S^{}M}}\sigma _0(B_{j_1}\mathrm{}B_{j_m})dxd\xi `$
$`+n^{d2}[{\displaystyle \frac{d1}{2}}(\left[M_m(\overline{j})\right]^2+(1+{\displaystyle \frac{\alpha }{2}})M_m(\overline{j})){\displaystyle _{S^{}M}}\sigma _0(B_{j_1}\mathrm{}B_{j_m})dxd\xi `$
$`+M_m(\overline{j}){\displaystyle _{S^{}M}}sub(B_{j_1}\mathrm{}B_{j_m})dxd\xi ]\}`$
$`+O(n^{d3}),n\mathrm{}.`$
Now the symmetrization argument from \[K, S, O1, GO2, GO3, O2\] comes into play. Note that the domain of summation over $`j_1,\mathrm{},j_m`$ in (1.12) is symmetric with respect to the permutations of the $`j`$’s. Also the factors involving

$$\sigma _0(B_{j_1}\mathrm{}B_{j_m})=\sigma _0(B_{j_1})\mathrm{}\sigma _0(B_{j_m})$$

are symmetric. Therefore one can write

$`{\displaystyle \underset{j_1+\mathrm{}+j_m=0}{}}`$ $`M_m(\overline{j}){\displaystyle _{S^{}M}}\sigma _0(B_{j_1}\mathrm{}B_{j_m})𝑑x𝑑\xi `$
$`={\displaystyle \underset{j_1+\mathrm{}+j_m=0}{}}\left({\displaystyle \frac{1}{m!}}{\displaystyle \underset{\tau S_m}{}}M_m(\overline{j}_\tau )\right){\displaystyle _{S^{}M}}\sigma _0(B_{j_1}\mathrm{}B_{j_m})𝑑x𝑑\xi ,`$
where $`S_m`$ is the set of all permutations of $`j_1,\mathrm{},j_m`$ and $`\overline{j}_\tau :=(j_{\tau _1},\mathrm{},j_{\tau _m})`$. After an application of the usual Hunt–Dyson combinatorial formula (see (3.4) in Section 3)

$$\underset{\tau S_m}{}M_m(\overline{j}_\tau )=\underset{\tau S_m}{}\underset{k=1}{\overset{m}{}}\frac{min(0,j_{\tau _1}+\mathrm{}+j_{\tau _k})}{k},$$

all summations but one in $`_{j_1+\mathrm{}+j_m=0}`$ become free and can be carried out.

Now we also see the two difficulties with the computation of the $`n^{d2}`$ term. The first difficulty is in computing the minimum raised to power $`2`$. Here one needs a generalization of the Hunt–Dyson formula, see Theorem 3.1. The second difficulty is that the factor

$$_{S^{}M}sub(B_{j_1}\mathrm{}B_{j_m})dxd\xi $$

is generally speaking not symmetric with respect to the permutations $`j_1,\mathrm{},j_m`$, and the straightforward symmetrization as above fails. Indeed, recall that

(1.13) $`sub(B_{j_1}\mathrm{}B_{j_m})`$ $`={\displaystyle \underset{k=1}{\overset{m}{}}}sub(B_{j_k}){\displaystyle \underset{\genfrac{}{}{0.0pt}{}{p=1}{pk}}{\overset{m}{}}}\sigma _0(B_{j_p})`$
$`+{\displaystyle \frac{1}{2i}}{\displaystyle \underset{1k<lm}{}}\{\sigma _0(B_{j_k}),\sigma _0(B_{j_l})\}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{p=1}{pk,pl}}{\overset{m}{}}}\sigma _0(B_{j_p}).`$
where $`\{,\}`$ stands for the Poisson bracket. We see that the first sum in (1.13) is symmetric, whereas the second one is generally speaking not (not even after the integration over $`S^{}M`$). However, each of the $`m(m1)/2`$ terms in the second sum on the right-hand side of (1.13) possesses a partial symmetry. More precisely, we are allowed to permute the principal symbols that do not enter the Poisson bracket. This fact can be utilized in a modification of the symmetrization procedure. The corresponding contribution is lengthy, this is where the functional $`\mathrm{{\rm Y}}`$ arises. ∎

###### Remark 1.6.

We would like to mention here that in the computation of $`\mathrm{{\rm Y}}`$ one really needs the information on the set of values  $`\left\{M_m(\overline{j}_\tau )\right\}_{\tau S_m}`$ counted with multiplicities, that is the Bohnenblust–Spitzer theorem (see Section 3 for its statement), and not just a formula for $`_{\tau S_m}[M_m(\overline{j}_\tau )]^p`$ for $`p=1,2`$.

## 2. Explicit asymptotic formulas for $`\mathrm{log}detP_nBP_n`$, as $`n\mathrm{}`$.

Theorem 1.1 and 1.2 give an expression for $`Tr\mathrm{log}P_nBP_n=\mathrm{log}detP_nBP_n`$ as a sum of $`TrP_n(\mathrm{log}B)P_n`$ and the two lower order corrections, as $`n\mathrm{}`$. We would like to compute the coefficients in the asymptotic expansion of $`\mathrm{log}detP_nBP_n`$, as $`n\mathrm{}`$.

An auxiliary asymptotic expansion for

(2.1)
$$TrP_nGP_n=\underset{k=1}{\overset{n}{}}Tr(\pi _kG),n\mathrm{},$$

where $`G\mathrm{\Psi }^0(M)`$, is given in Proposition 2.1 below. A corresponding result for the coefficients of $`n^d`$, $`n^{d1}`$ and $`\mathrm{log}n`$ for $`d2`$ can be found in \[GO2, after Lemma 0.2\]. To prove Proposition 2.1, we sum over $`k=1,\mathrm{},n`$ in (1.11), see \[Gi2\] for details. The subtle point is the constant coefficient in (2.1), which we need for dimension $`d=1,2`$. The terms of all orders in (1.11), and also the possible rapidly decaying term, will contribute to it. Let us therefore for $`d=1,2`$ make an additional assumption

(2.2)
$$\underset{l=0}{\overset{\mathrm{}}{}}\left|R_l(G)\right|<\mathrm{}$$

and set, for all $`k`$,

(2.3)
$$ϵ_k(G):=Tr(\pi _kG)\underset{l=0}{\overset{+\mathrm{}}{}}k^{d1l}R_l(G).$$

If (2.2) holds, the series in (2.3) (and also in (2.4) below) is absolutely convergent, and also for any $`N`$ there exists $`c_N<\mathrm{}`$ such that $`|ϵ_k|c_Nk^N`$, $`k`$. Set

(2.4)
$$C(G):=\underset{k=1}{\overset{\mathrm{}}{}}ϵ_k(G).$$

Note that in (2.4) and also in the proposition below there only appears the series $`_{l=0}^{\mathrm{}}R_l(G)`$. However we need the absolute convergence (2.2) in the proof of the remainder estimate.

Let $`\gamma `$ denote the Euler constant and $`\zeta `$ the Riemann zeta function.

###### Proposition 2.1.

Let $`M`$ be a Zoll manifold of dimension $`d`$. Let $`P_n`$ be as above and assume that $`G\mathrm{\Psi }^0(M)`$. For $`d=1,2`$, assume in addition that (2.2) holds, and let $`C(G)`$ be defined by (2.4). Then the following holds, as $`n\mathrm{}`$,
(i) for $`d=1`$,

$`Tr`$ $`P_nGP_n=nR_0(G)+\mathrm{log}nR_1(G)`$
$`+\left(C(G)+\gamma R_1(G)+{\displaystyle \underset{l=2}{\overset{\mathrm{}}{}}}\zeta (l)R_l(G)\right)+{\displaystyle \frac{1}{n}}\left({\displaystyle \frac{1}{2}}R_1(G)R_2(G)\right)+O\left({\displaystyle \frac{1}{n^2}}\right),`$
(ii) for $`d=2`$,

$`TrP_nGP_n`$ $`=n^2{\displaystyle \frac{1}{2}}R_0(G)+n\left({\displaystyle \frac{1}{2}}R_0(G)+R_1(G)\right)+\mathrm{log}nR_2(G)`$
$`+\left(C(G)+\gamma R_2(G)+{\displaystyle \underset{l=2}{\overset{\mathrm{}}{}}}\zeta (l)R_{l+1}(G)\right)+O\left({\displaystyle \frac{1}{n}}\right),`$
(iii) for $`d3`$,

$`Tr`$ $`P_nGP_n=n^d{\displaystyle \frac{1}{d}}R_0(G)+n^{d1}\left({\displaystyle \frac{1}{2}}R_0(G)+{\displaystyle \frac{1}{d1}}R_1(G)\right)`$
$`+n^{d2}\left({\displaystyle \frac{d1}{12}}R_0(G)+{\displaystyle \frac{1}{2}}R_1(G)+{\displaystyle \frac{1}{d2}}R_2(G)\right)+\mathrm{log}nR_d(G)+O(n^{d3}).`$
###### Remark 2.1.

We see that $`TrP_n(\mathrm{log}B)P_n`$ in Theorem 1.1 and 1.2 ($`G=\mathrm{log}B`$) contributes to the leading asymptotic term of order $`n^d`$, and also to all lower order terms of order $`n^j`$, $`j=d1,\mathrm{},1,0,1,\mathrm{}`$, and to the logarithmic term $`\mathrm{log}n`$, as $`n\mathrm{}`$. In the classical SSLT the situation is much simpler: $`\mathrm{log}B`$ is just the Toeplitz matrix of the operator of multiplication by $`\mathrm{log}b`$, and so $`TrP_n(\mathrm{log}B)P_n=(2n+1)\widehat{(\mathrm{log}b)}_0`$.

Now we are ready to state the two corollaries.

###### Corollary 2.2.

Let $`B\mathrm{\Psi }^0(𝕊^1)`$ have a strictly positive principal symbol and be such that a certain symbolic norm of $`IB`$ is sufficiently small. Assume also that (2.2) holds. Then the following holds, as $`n\mathrm{}`$,

$$\mathrm{log}detP_nBP_n=c_1n+c_{\mathrm{log}}\mathrm{log}n+c_0+c_1\frac{1}{n}+O\left(\frac{1}{n^2}\right),$$

where the coefficients are the sums of the corresponding coefficients from Theorem 1.1 and Proposition 2.1(i).

Assume further that $`\sigma _0(B)`$ and $`sub(B)`$ do not depend on the direction of $`\xi `$, that is $`\sigma _0(B)(x,\xi )=b_0(x)`$ and $`sub(B)(x,\xi )=b_{sub}(x)|\xi |^1`$, for $`(x,\xi )S^{}𝕊^1`$. Assume also that $`b_2=0`$. Then the following holds, as $`n\mathrm{}`$,

(2.5) $`\mathrm{log}`$ $`detP_nBP_n=n2{\displaystyle _0^{2\pi }}\mathrm{log}b_0(x){\displaystyle \frac{dx}{2\pi }}`$
$`+\mathrm{log}n2{\displaystyle _0^{2\pi }}{\displaystyle \frac{b_{\mathrm{sub}}(x)}{b_0(x)}}{\displaystyle \frac{dx}{2\pi }}`$
$`+\left({\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\widehat{(\mathrm{log}b_0)}_k+C(\mathrm{log}B)+\gamma R_1(\mathrm{log}B)+{\displaystyle \underset{l=2}{\overset{\mathrm{}}{}}}\zeta (l)R_l(\mathrm{log}B)\right)`$
$`+{\displaystyle \frac{1}{n}}\left({\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k\widehat{(\mathrm{log}b_0)}_k\left(\widehat{b_{\mathrm{sub}}/b_0}\right)_k+{\displaystyle _0^{2\pi }}\left[{\displaystyle \frac{b_{\mathrm{sub}}(x)}{b_0(x)}}+\left({\displaystyle \frac{b_{\mathrm{sub}}(x)}{b_0(x)}}\right)^2\right]{\displaystyle \frac{dx}{2\pi }}\right)`$
$`+O\left({\displaystyle \frac{1}{n^2}}\right),`$
where $`C(\mathrm{log}B)`$ is given by (2.4).

###### Remark 2.2.

In some simple cases, for instance for $`b_{\mathrm{sub}}(x)=\pm \frac{1}{2}b_0(x)`$, the left-hand side in (2.5) can be computed explicitly. The coefficients of $`n`$, $`\mathrm{log}n`$, and $`\frac{1}{n}`$ on the right in (2.5) in these cases are as expected, see also Remark 2.7.

###### Corollary 2.3.

Let $`M`$ be a Zoll manifold of dimension $`d2`$. Assume that $`P_n`$ and $`A`$ are as in Theorem 1.2. Let $`B\mathrm{\Psi }^0(M)`$ have a strictly positive principal symbol and be such that the symbolic norm of $`IB`$ is sufficiently small. For $`d=2`$, assume in addition (2.2). Then the following holds, as $`n\mathrm{}`$,

(2.6)
$$\mathrm{log}detP_nBP_n=C_d^{(d)}n^d+C_{d1}^{(d)}n^{d1}+C_{d2}^{(d)}n^{d2}+C_{\mathrm{log}}^{(d)}\mathrm{log}n+O\left(n^{d3}\right),$$

where the coefficients are the sums of the corresponding coefficients from Theorem 1.2 and Proposition 2.1(ii) or (iii). If one counts the logarithmic term, this expansion is fourth order for $`d=2,3`$ and third order for $`d4`$.

###### Remark 2.3.

The coefficients $`C_d^{(d)}`$, $`C_{d1}^{(d)}`$, and also $`C_{\mathrm{log}}^{(d)}`$, $`d`$, have been found in \[GO1, GO2\].

###### Remark 2.4.

The most interesting coefficient in (2.6) is the constant one, since one can think of $`\mathrm{exp}C_0^{(d)}`$ as of a regularized determinant of $`B\mathrm{\Psi }^0(M)`$, see \[GO1, O2\]. The sum

$$\gamma R_d(\mathrm{log}B)+\underset{l=2}{\overset{\mathrm{}}{}}\zeta (l)R_{l+d1}(\mathrm{log}B)$$

will for all $`d`$ be a part of $`C_0^{(d)}`$. For $`d=1`$, Corollary 2.2 gives a full expression for $`C_0^{(1)}`$. For $`d=2`$, Corollary 2.3 gives a full expression for $`C_0^{(2)}`$, which is quite lengthy.

###### Remark 2.5.

Let us compare the result of Corollary 2.2 with a generalization of SSLT to the case of $`B`$ being an operator of multiplication by a function $`b(x)`$ having discontinuities which is due to H. Widom and E. Basor. In this case $`\mathrm{log}b(x)`$ also has discontinuities, and so the series $`_k|k||\widehat{(\mathrm{log}b)}_k|^2`$ diverges logarithmically. The following third order asymptotic formula holds for the operator of multiplication by a piecewise $`C^2`$ function $`b(x)`$

(2.7)
$$\mathrm{log}detP_nBP_n=a_1n+a_2\mathrm{log}n+a_3+o(1),n\mathrm{},$$

where $`a_1`$ as in (2.5), the coefficient $`a_2`$ has been computed by H. Widom in \[W2\], and the constant term $`a_3`$ has been found by E. Basor in \[B\]. Note that the matrix of $`B`$ in (2.7) is still Toeplitz, the logarithmic order of the subleading term being due to a slower decay of the Fourier coefficients of $`b(x)`$. In our case the matrix of the operator $`B\mathrm{\Psi }^0(𝕊^1)`$ is not Toeplitz (see Remark 2.7), and the logarithm comes from the contribution of $`sub(B)`$.

###### Remark 2.6.

It would be interesting to find a compact formula for the constant term in (2.5). We mention that the constant $`a_3`$ in (2.7) found in \[B\] has a form similar to the one in (2.5). It contains a “finite” term and and an infinite series of certain integrals multiplied by the values of the Riemann zeta function at the points $`3,5,\mathrm{}`$. Interestingly, an “invariant” form of that series has been found in \[W4\]. It is written as a single integral involving the function

$$\mathrm{\Psi }(x):=\frac{d}{dx}\mathrm{log}\mathrm{\Gamma }(x).$$

This gives the hope that a similar formula can be found for the constant (2.5). We have done some computations trying to find the constant term in (2.5), and the function $`\mathrm{\Psi }(x)`$ has been appearing there.

###### Remark 2.7.

The matrix interpretation of Corollary 2.2 is as follows. Assume for simplicity that $`B\mathrm{\Psi }^0(𝕊^1)`$ is as in the second part of Corollary 2.2, that is $`\sigma _0(B)(x,\xi )=b_0(x)`$ and $`sub(B)(x,\xi )=b_{sub}(x)|\xi |^1`$, for all $`(x,\xi )S^{}𝕊^1`$. Assume also that $`b_2=b_3=\mathrm{}=0`$. Let $`B_0`$ and $`B_{\mathrm{sub}}`$ be the operators of multiplication by $`b_0`$ and $`b_{\mathrm{sub}}`$, respectively. Let $`D`$ be the linear operator in $`L^2(𝕊^1)`$ such that

$$De^{ikx}=\{\begin{array}{cc}\frac{1}{|k|}e^{ikx},\hfill & |k|1\hfill \\ 0,\hfill & k=0.\hfill \end{array}$$

Note that this is not a differential, but rather a smoothing operator of order $`1`$. There is known a correspondence between the classical PsDO’s on the circle and their discrete counterparts, see \[M\] for details. By that correspondence, the zeroth order PsDO $`B`$ we started with, equals $`B_0+B_{\mathrm{sub}}D`$. Introduce two Toeplitz matrices, $`\widehat{B}_0:=\{\widehat{(b_0)}_{jk}\}_{j,k}`$ and $`\widehat{B}_{\mathrm{sub}}:=\{\widehat{(b_{\mathrm{sub}})}_{jk}\}_{j,k}`$. Set also $`\widehat{D}:=diag(\mathrm{},\frac{1}{3},\frac{1}{2},1,0,1,\frac{1}{2},\frac{1}{3},\mathrm{})`$. Then the matrix representation of $`B_0+B_{\mathrm{sub}}D`$ is $`\widehat{B}_0+\widehat{B}_{\mathrm{sub}}\widehat{D}`$. Finally, set $`\widehat{P}_n=diag(\mathrm{},0,1,\mathrm{},1,0,\mathrm{})`$ ($`2n+1`$ ones). We see that Corollary 2.2 gives a fourth order asymptotics of the determinant of the truncated matrix $`\widehat{P}_n(\widehat{B}_0+\widehat{B}_{\mathrm{sub}}\widehat{D})\widehat{P}_n`$.

Now we can reformulate the question of finding the constant term in (2.5) in purely matrix terms. Drop the hats and the dots for brevity. Let $`C_1`$ be a Toeplitz matrix that corresponds to the operator of multiplication by $`b_{\mathrm{sub}}/b_0`$, and let the matrix $`D`$ be as above. Clearly, the matrices $`C_1`$ and $`D`$ do not commute. Assume that the matrix $`\mathrm{log}(IC_1D)`$ is well-defined. The question is to compute the constant coefficient in $`TrP_n\mathrm{log}(IC_1D)P_n`$, or which is the same, the constant coefficient in

(2.8)
$$TrP_n\mathrm{log}(ID^{1/2}C_1D^{1/2})P_n,n\mathrm{}.$$

As we have noticed in Remark 2.1, this question is trivial for a Toeplitz matrix $`T`$ in place of $`D^{1/2}C_1D^{1/2}`$.

## 3. Generalized Hunt–Dyson combinatorial formula

In this section we state the generalized Hunt–Dyson formula (gHD) and remind the reader the Bohnenblust–Spitzer theorem (BSt). After that we briefly explain how the gHD is derived from the BSt in \[Gi3\]. We refer to \[Gi3, Section 4\] for the first step of an independent proof of the gHD. This first step is a generalization of F. J. Dyson’s idea on which the proof of the usual Hunt–Dyson formula in \[K\] is based. It is also explained in \[Gi3, Section 4\] how to reprove the BSt starting with the gHD. See \[Gi1, Chapter 2\] for the details of the independent proof of the gHD.

We state the result for the maximum, for the corresponding result for the minimum one should replace the positive parts with the negative parts. For $`a`$, $`n`$ denote

$$(a)_+:=\mathrm{max}(0,a),(a)_+^n:=((a)_+)^n.$$

Fix any $`m`$ and assume $`a_1,\mathrm{},a_m`$. Let $`S_m`$ be the set of all permutations $`\tau `$ of the numbers $`1,\mathrm{},m`$. For each $`\tau =\left(\genfrac{}{}{0.0pt}{}{1,2,\mathrm{},m}{\tau _1,\tau _2,\mathrm{},\tau _m}\right)S_m`$ denote

$$a_\tau :=(a_{\tau _1},\mathrm{},a_{\tau _m}).$$

Introduce the notation

(3.1)
$$M_j(a_\tau ):=\{\begin{array}{cc}\mathrm{max}(0,a_{\tau _1},a_{\tau _1}+a_{\tau _2},\mathrm{},a_{\tau _1}+\mathrm{}+a_{\tau _j}),\hfill & j=1,\mathrm{},m,\hfill \\ 0,\hfill & j=0.\hfill \end{array}$$

Fix any $`j=1,\mathrm{},m`$. For arbitrary

$$k_11,\mathrm{},k_j1,k_1+\mathrm{}+k_j=m,$$

we introduce the notation

(3.2) $`k_1(a_\tau )`$ $`:=a_{\tau _1}+\mathrm{}+a_{\tau _{k_1}}`$
$`k_2(a_\tau )`$ $`:=a_{\tau _{k_1+1}}+\mathrm{}+a_{\tau _{k_1+k_2}}`$
$`\mathrm{}`$
$`k_j(a_\tau )`$ $`:=a_{\tau _{k_1+\mathrm{}+k_{j1}+1}}+\mathrm{}+a_{\tau _{k_1+\mathrm{}+k_{j1}+k_j}}.`$
Each of $`k_l(a_\tau )`$, $`l=1,\mathrm{},j`$, is a sum of $`k_l`$ permuted variables out of $`a_{\tau _1},\mathrm{},a_{\tau _m}`$ so that each of the permuted variables enters exactly one sum. Recall that $`\left(\genfrac{}{}{0.0pt}{}{n}{l_1,\mathrm{},l_j}\right):=\frac{n!}{l_1!\mathrm{}l_j!}`$ denotes a multinomial coefficient, here $`n,j`$, $`l_1,\mathrm{},l_j\{0\}`$, and $`l_1+\mathrm{}+l_j=n`$. We are ready to state the generalized Hunt–Dyson formula (gHD).

###### Theorem 3.1.

For an arbitrary power $`n`$, an arbitrary number of variables $`m`$, and for arbitrary $`a_1,\mathrm{},a_m`$, the following holds

(3.3) $`{\displaystyle \underset{\tau S_m}{}}`$ $`\left[(M_m(a_\tau ))^n(M_{m1}(a_\tau ))^n\right]`$
$`={\displaystyle \underset{\tau S_m}{}}{\displaystyle \underset{j=1}{\overset{min(m,n)}{}}}{\displaystyle \frac{1}{j!}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{k_1,\mathrm{},k_j1}{k_1+\mathrm{}+k_j=m}}{}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{l_1,\mathrm{},l_j1}{l_1+\mathrm{}+l_j=n}}{}}\left({\displaystyle \genfrac{}{}{0.0pt}{}{n}{l_1,\mathrm{},l_j}}\right){\displaystyle \frac{(k_1(a_\tau ))_+^{l_1}}{k_1}}\mathrm{}{\displaystyle \frac{(k_j(a_\tau ))_+^{l_j}}{k_j}}.`$
###### Remark 3.1.

In the case $`n=1`$ we obtain $`j=1,k_1=m,l_1=1,`$ and (3.3) becomes the usual Hunt–Dyson combinatorial formula (HD) \[K, after (4.8)\]

(3.4) $`{\displaystyle \underset{\tau S_m}{}}\left[M_m(a_\tau )M_{m1}(a_\tau )\right]`$ $`={\displaystyle \underset{\tau S_m}{}}{\displaystyle \frac{(a_{\tau _1}+\mathrm{}+a_{\tau _m})_+}{m}}`$
$`=(m1)!(a_1+\mathrm{}+a_m)_+.`$
Recall now the statement of the Bohnenblust–Spitzer theorem (BSt) \[S, Theorem 2.2\]. It asserts that for any $`m`$ and arbitrary $`a_1,\mathrm{},a_m`$, the set $`\{M_m(a_\tau )\}_{\tau S_m}`$ contains the same numbers with the same multiplicities as the set of sums of positive parts of the sums of $`a_1,\mathrm{},a_m`$, arranged according to the cyclic representations of all $`m!`$ permutations. This becomes clear if we consider a simple example. Let us choose $`m=3`$ and any $`a_1,a_2,a_3`$. The symmetric group $`S_3`$ consists of six permutations that can be written via the cyclic representations as

$$S_3=\{(123),(132),(12)(3),(13)(2),(23)(1),(1)(2)(3)\}$$

In this case the BSt states that the set

(3.5)
$$\left\{\mathrm{max}(0,a_{\sigma _1},a_{\sigma _1}+a_{\sigma _2},a_{\sigma _1}+a_{\sigma _2}+a_{\sigma _3})\right\}_{\sigma S_3}$$

contains the same numbers with the same multiplicities as the set

(3.6) $`\{`$ $`(a_1+a_2+a_3)_+,(a_1+a_3+a_2)_+,`$
$`(a_1+a_2)_++(a_3)_+,(a_1+a_3)_++(a_2)_+,(a_2+a_3)_++(a_1)_+,`$
$`(a_1)_++(a_2)_++(a_3)_+\}.`$
Note that a certain maximum of zero and accumulating sums of the permuted variables does not need to equal the element of the set on the right-hand side corresponding to the cyclic representation of that permutation. The statement of the BSt is merely that the whole multisets are identical.

###### Proof of Theorem 3.1: an outline.

By the BSt, for any $`m`$ and arbitrary real $`a_1,\mathrm{},a_m`$, the counterparts of the sets (3.5) and (3.6) contain the same numbers with the same multiplicities. Therefore for any function $`f:_+`$, the sum $`_{\tau S_m}f(M_m(a_\sigma ))`$ equals the sum of $`f`$’s values over the counterpart of (3.6).

It turns out that for the multinomial function $`f(t)=t^m`$ for any $`n`$, a further computation can be performed which leads to the gHD, see \[Gi3, Section 2\].

Without going into details let us just mention here that when we apply the multinomial formula in the sum over the analog of (3.6), some terms have zero powers. Recall that we want all factors to be present in the right-hand side of the gHD (3.3). However it turns out that the terms having at least one zero power can be summed together, and their sum gives exactly the sum of the $`n`$th power of the “previous” maximum $`M_{m1}(a_\sigma )`$ over $`S_m`$.

At this step we use the principle of inclusion and exclusion and the Cauchy and Cayley identities from the theory of partitions. ∎

###### Remark 3.2.

The steps of the derivation of the gHD from the BSt can be reversed. After that having started with the gHD, we can conclude that for any monomial $`f(t)`$ its sum over (3.5) equals the sum over (3.6). Due to an additional linearity, this actually holds for an arbitrary polynomial $`f(t)`$. Now using the polynomial interpolation we arrive at the BSt, see \[Gi3, Section 4\].

###### Remark 3.3.

The independent proof of the gHD in \[Gi1, Chapter 2\] proceeds by induction on the power $`n`$. The base of induction is the usual HD ($`n=1`$). In the proof of the inductive step, the key cancellation of the highest power $`n`$ of the maximum after taking a sum over $`S_m`$ follows from a generalization of F. J. Dyson’s argument from the proof of the usual Hunt–Dyson formula in \[K\], see \[Gi3, Section 4\]. The proof of the inductive step is however quite technical.

###### Remark 3.4.

In \[RS\], the authors rediscover a version of the usual Hunt–Dyson formula starting from the much more powerful BSt.

Acknowledgments. Theorem 1.2 and Theorem 3.1 appear in the author’s Ph.D. thesis. I would like to express a deep gratitude to my thesis adviser Ari Laptev, for suggesting the problem and his constant attention to the work.

I would also like to thank Percy Deift for giving a reference \[RS\] from which I learned about \[S\] and the Bohnenblust–Spitzer combinatorial theorem, shortly after having discovered and proved the generalized Hunt–Dyson formula.

Finally, I would like to thank the Department of Mathematics of the Royal Institute of Technology (KTH), Stockholm, for excellent working conditions and a generous financial support during the whole period of my graduate studies.


text:
# Momentum dependence of symmetry potential in asymmetric nuclear matter for transport model calculations

## I Introduction

Momentum dependent mean-fields for transport model calculations of heavy ion collisions have been in usage for many years now . So far the attention has been for a momentum dependent potential which does not distinguish between neutrons and protons. This is adequate for systems which have $`NZ`$ where $`N`$ is the number of neutrons and $`Z`$ the number of protons. One main focus of nuclear physics research today is to explore both structures of exotic nuclei in regions far-off the stability line and novel properties of neutron-rich nuclear matter. The latter can be investigated by using collisions induced by neutron-rich nuclei at intermediate to high energies. To interpret critically data from these collisions and to extract accurately properties of neutron-rich nuclear matter, advanced transport model calculations are necessary. In asymmetric nuclear matter, the one body potential seen by a proton is different from that seen by a neutron. This has been implemented in BUU (Boltzmann-Uehling-Uhlenbeck) calculations but with a simplification that the potentials are taken to be momentum independent . The present work aims to correct this deficiency. That the momentum dependence will be different for neutrons and protons is of course well-known and has been the subject of quite sophisticated many body calculations, see, e.g, ref.  for a recent review. We do not aim to add anything fundamental in this regard. Our objective is to obtain a parametrized version which displays the main characteristics of momentum dependence in asymmetric matter and is still usable in practical BUU calculations. Major advances in this direction were already made: see articles by Bombaci  and Prakash et al.. We add to this. We will not only extend the simplest momentum dependent potential  to include isospin but also extend the improved treatment  subsequently introduced to include isospin. Thus this is an extension of the work reported in .

## II A momentum dependent potential from a phenomenological interaction

An effective momentum dependent potential can be deduced from phenomenological interactions. We take the Gogny interaction  to obtain an idea of the momentum dependence. There are many reasons for this choice. It has been used in detailed fits for spectra in finite nuclei. It gives accepted values for binding energy, saturation density, compressibility and symmetry energy in nuclear matter. It has been verified already  that the interaction produces a reasonable parameterisation for the real part of the optical potential in nuclear matter as a function of incident energy. The Skyrme interaction has a wrong asymptotic behaviour as a function of energy. (This is amplified in .) Since we want to devise a momentum dependence which should hold for beam energy as high as 1 GeV/nucleon (this would allow investigation of symmetry energy at higher than normal nuclear density) we discard the Skyrme interaction.

For the purpose of this work we will define nuclear matter to be an infinite system but without the restriction $`N=Z`$. Using the Gogny interaction, we deduce $`U(\rho ,\delta ,p,\tau )`$, the one body potential a particle of momentum $`p`$ and isospin $`\tau `$ feels in cold nuclear matter with density $`\rho `$ and asymmetry $`\delta \frac{\rho _n\rho _p}{\rho _n+\rho _p}`$. One then generalises to $`U`$ in the case of heavy ion collisions. For BUU calculations $`U`$ is the only quantity needed. But it is useful to also have an expression for $`V(\rho ,\delta )`$, the potential energy density in cold nuclear matter with a given density $`\rho `$ and asymmetry $`\delta `$. This allows one to deduce $`E/A`$ as a function of $`\rho `$ and $`\delta `$ which is, of course, of importance. The expression for $`V(\rho ,\delta )`$ can also be generalised to the case of heavy ion collisions and can be used to check, for example, the accuracy of energy conservation in a BUU simulation.

Fig.1: Single particle potential $`U`$ with respect to $`k`$ and total single particle energy ($`e`$) for neutron and proton for $`\delta `$=0.4. The solid line is the single particle potential for symmetric matter; here $`\rho =0.16fm^3`$.

Normally nuclear matter denotes an infinite nucleus with $`N=Z`$. Total potential energy in cold matter is deduced from

$`V_T`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{p_1\sigma _1\tau _1p_2\sigma _2\tau _2}{}}\stackrel{}{p}_1,\sigma _1,\tau _1,\stackrel{}{p}_2,\sigma _2,\tau _2|v(r)`$ (2)
$`\left(|\stackrel{}{p}_1,\sigma _1,\tau _1,\stackrel{}{p}_2,\sigma _2,\tau _2|\stackrel{}{p}_2,\sigma _2,\tau _2,\stackrel{}{p}_1,\sigma _1,\tau _1\right)`$
where,

$`v(r)`$ $`=`$ $`{\displaystyle \underset{i=1,2}{}}(W+BP_\sigma HP_\tau MP_\sigma P_\tau )_ie^{r^2/\mu _i^2}`$ (4)
$`+t_0(1+P_\sigma )\rho ^\alpha ({\displaystyle \frac{\stackrel{}{r}_1+\stackrel{}{r}_2}{2}})\delta (\stackrel{}{r}_1\stackrel{}{r}_2)`$
There are two finite range Gaussians and a density dependent zero-range force. The values of the parameters are given in . The one body potential $`U(\rho ,\delta ,p,\tau )`$ is obtained from

$`U(\rho ,\delta ,p,\tau )=_{p^{}\sigma ^{}\tau ^{}}\stackrel{}{p},\sigma ,\tau ,\stackrel{}{p}{}_{}{}^{},\sigma ^{},\tau ^{}|v(r)(|\stackrel{}{p},\sigma ,\tau ,\stackrel{}{p}{}_{}{}^{},\sigma ^{},\tau ^{}|\stackrel{}{p}{}_{}{}^{},\sigma ^{},\tau ^{},\stackrel{}{p},\sigma ,\tau )`$ plus rearrangement term which for nuclear matter is $`(3/2)t_0\alpha \rho ^{\alpha 1}(1/4)\rho ^2(1\delta ^2)`$.

The momentum dependence in $`U`$ comes entirely from the exchange term of the finite range part, i.e., from $`\stackrel{}{p},\stackrel{}{p}{}_{}{}^{}|e^{r^2/\mu ^2}|\stackrel{}{p}{}_{}{}^{},\stackrel{}{p}`$. Except for the momentum dependent part, very simple expressions for $`U`$ and $`V_T/A`$, the potential energy per particle are obtained for the Gogny potential. Thus $`U(\rho ,\delta ,p,\tau )=X+Y+Z`$ where $`X`$ arises from the direct term of the finite range interaction, $`Y`$ arises from the $`t_0`$ term (density dependent two body term) and $`Z`$ from the exchange term of the finite range interactions. For a given $`\tau `$, these are:

$`X`$ $`=`$ $`\rho \left({\displaystyle \underset{i=1,2}{}}\pi ^{3/2}\mu _i^3(W+B/2)_i\right)\rho _\tau \left({\displaystyle \underset{i=1,2}{}}\pi ^{3/2}\mu _i^3(H+M/2)_i\right)`$ (5)
$`Y`$ $`=`$ $`{\displaystyle \frac{3}{2}}t_0\rho ^\alpha (\rho \rho _\tau )+{\displaystyle \frac{3}{2}}t_0\alpha \rho ^{\alpha 1}{\displaystyle \frac{1}{4}}\rho ^2(1\delta ^2)`$ (6)
$`\mathrm{and}Z`$ $`=`$ $`\sqrt{\pi }\left[{\displaystyle \underset{i=1,2}{}}Z_i(p,\tau )(W2B+H+2M)_i+{\displaystyle \underset{i}{}}Z_i(p,\tau ^{})(H+2M)_i\right].`$ (7)
where,

$`Z_i(p,\tau )`$ $`=`$ $`{\displaystyle \frac{1}{\mu _ik}}\left[e^{(\mu _i(k_F(\tau )k)/2)^2}e^{(\mu _i(k_F(\tau )k)/2)^2}\right]`$ (9)
$`+{\displaystyle \frac{\sqrt{\pi }}{2}}[erf({\displaystyle \frac{\mu _i}{2}}(k_F(\tau )k)+erf({\displaystyle \frac{\mu _i}{2}}(k_F(\tau )+k)].`$
Here $`\tau ^{}\tau `$, the isospin of the particle whose one body potential is being sought.

Similarly, $`V_T/A`$=potential energy per particle, has contributions from the direct term of the finite range force, from the density dependent $`t_0`$ term and the exchange term of the finite range force. Denoting the first two by $`X^{}`$ and $`Y^{}`$ respectively, explicit expressions for these are:

$`X^{}`$ $`=`$ $`\rho [{\displaystyle \underset{i=1,2}{}}\pi ^{3/2}\mu _i^3({\displaystyle \frac{W}{2}}+{\displaystyle \frac{B}{4}}{\displaystyle \frac{H}{4}}{\displaystyle \frac{M}{8}})_i]\rho \delta ^2[{\displaystyle \underset{i=1,2}{}}({\displaystyle \frac{H}{4}}+{\displaystyle \frac{M}{8}})_i]`$ (10)
$`Y^{}`$ $`=`$ $`{\displaystyle \frac{3}{8}}t_0\rho ^{\alpha +1}(1\delta ^2).`$ (11)
Fig.2 $`E_{sym}(\rho )`$ as a function of $`\rho `$.

We do not write down the explicit expression for $`Z^{}`$. It is obtained from the expression of $`Z`$ above after a further integration over $`p`$, sum over $`\sigma ,\tau `$ and dividing the answer by 2.

The one body potentials as a function of $`k`$ for neutrons and protons in cold matter as predicted by the Gogny potential for $`\delta =0`$ and $`\delta =0.4`$ are shown in Fig. 1. They are quite similar to other calculations of $`U(\rho ,\delta ,p,\tau )`$ available in the literature. We compare, in particular, to the $`U(\rho ,\delta ,p,\tau )`$ given in Fig.4 in . That figure is for $`\rho =.17`$ and obtained from Brueckner-Hartee-Fock calculation. With Gogny interaction we find that the equation of state (EOS) of asymmetric nuclear matter can be written as $`\frac{E}{A}(\rho ,\delta )\frac{E}{A}(\rho ,0)+E_{sym}(\rho )\delta ^2`$, in agreement with the empirical parabolic law found by all many body theories. Different calculations depending on the many body approaches and the interactions used, give very different behaviors for $`E_{sym}(\rho )`$, especially at high densities. In some, such as, the Relativistic Mean Field model and the Brueckner-Hartree-Fock, it is a continuously rising function of $`\rho `$, In others, such as, the Variational Many-Body Approach, it rises in the beginning and then begins to fall. Within the Hartree-Fock approach using all 86 Skyrme effective interactions widely used currently in nuclear structure studies, it was found that about $`1/3(2/3)`$ of them lead to symmetry energies in the first (second) category. Gogny interaction with the default parameters has behaviour in the second category. It is interesting to see what causes the fall and how the parameters of the interaction can be changed minimally to alter this behaviour. In Fig.2 we plot contributions to $`E_{sym}(\rho )`$ from the Gogny interaction and we show separately the term coming from (a) the direct term of the finite range part, (b) the exchange term of the finite range part and (c) the density dependent zero range part. It is the latter that causes the bending over. For example, if we choose the density dependent term to be $`\delta `$ independent, $`E_{sym}(\rho )`$ will continue to rise.

Lastly, we write the potential energy density due to Gogny interaction in a form that is common practice in BUU literature. Thus

$`V(\rho ,\delta )`$ $`=`$ $`52.41({\displaystyle \frac{\rho ^2}{\rho _0}})+37.70({\displaystyle \frac{\rho ^2}{\rho _0}})\delta ^2+{\displaystyle \frac{102.6}{\sigma +1}}{\displaystyle \frac{\rho ^{\sigma +1}}{\rho _0^\sigma }}(1\delta ^2)`$ (14)
$`+{\displaystyle \frac{1}{\rho _0}}{\displaystyle \underset{\tau ,\tau ^{}}{}}{\displaystyle }{\displaystyle }f_\tau (\stackrel{}{r},\stackrel{}{p})f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})[(41.51+46.02\delta _{\tau ,\tau ^{}})e^{(\stackrel{}{k}\stackrel{}{k}^{})^2\mu _1^2/4}`$
$`+(38.62+17.25\delta _{\tau \tau ^{}})e^{(\stackrel{}{k}\stackrel{}{k}^{})^2\mu _2^2/4}]d^3pd^3p^{}.`$
Here, as in previous work, all quoted numbers are in MeV, $`\rho _0=0.16fm^3`$. Also $`\sigma =\alpha `$(Gogny)+1=4/3.

Given that the momentum dependence generated by the Gogny potential comes from two Gaussians, we proceed to find a simpler version for momentum dependent potential for asymmetric matter. We will extend the parametrisation of . Subsequently we will find an even simpler version, the kind that was used in first applications of momentum dependent potentials for heavy ion collisions.

## III A simple momentum dependent potential for BUU calculations

The simplest generalisation of the potential energy density of Eq.(5.4) of  to asymmetric nuclear matter is

$`V(\rho ,\delta )`$ $`=`$ $`{\displaystyle \frac{A_1}{2\rho _0}}\rho ^2+{\displaystyle \frac{A_2}{2\rho _0}}\rho ^2\delta ^2+{\displaystyle \frac{B}{\sigma +1}}{\displaystyle \frac{\rho ^{\sigma +1}}{\rho _0^\sigma }}(1x\delta ^2)`$ (16)
$`+{\displaystyle \frac{1}{\rho _0}}{\displaystyle \underset{\tau ,\tau ^{}}{}}C_{\tau ,\tau ^{}}{\displaystyle d^3pd^3p^{}\frac{f_\tau (\stackrel{}{r},\stackrel{}{p})f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}`$
The parameter $`x`$ is introduced to cover the largely uncertain behavior of nuclear symmetry energy $`E_{sym}(\rho )`$ as discussed in the previous section. For the choice of $`x=1`$ (same as in Gogny) in the term containig $`B`$, the symmetry energy will bend over beyond a density $`\rho `$; for the choice $`x=0`$ the symmetry energy will continue to rise with density. In the above, $`C_{1/2,1/2}=C_{1/2,1/2}=C_{like}`$ and $`C_{1/2,1/2}=C_{1/2,1/2}=C_{unlike}`$. In terms of interactions between like and unlike particles, the above equation is equivalent to

$`V(\rho ,\delta )`$ $`=`$ $`{\displaystyle \frac{A_u\rho _n\rho _p}{\rho _0}}+{\displaystyle \frac{A_l}{2\rho _0}}(\rho _n^2+\rho _p^2)+{\displaystyle \frac{B}{\sigma +1}}{\displaystyle \frac{\rho ^{\sigma +1}}{\rho _0^\sigma }}(1x\delta ^2)`$ (18)
$`+{\displaystyle \frac{1}{\rho _0}}{\displaystyle \underset{\tau ,\tau ^{}}{}}C_{\tau ,\tau ^{}}{\displaystyle d^3pd^3p^{}\frac{f_\tau (\stackrel{}{r},\stackrel{}{p})f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}`$
where $`A_1=(A_u+A_l)/2`$ and $`A_2=(A_lA_u)/2`$. The one-body potential needed for BUU computations is given by

$`U(\rho ,\delta ,\stackrel{}{p},\tau )`$ $`=`$ $`A_u{\displaystyle \frac{\rho _\tau ^{}}{\rho _0}}+A_l{\displaystyle \frac{\rho _\tau }{\rho _0}}+B({\displaystyle \frac{\rho }{\rho _0}})^\sigma (1x\delta ^2)x{\displaystyle \frac{B}{\sigma +1}}{\displaystyle \frac{\rho ^{\sigma +1}}{\rho _0^\sigma }}{\displaystyle \frac{d\delta ^2}{d\rho _\tau }}`$ (20)
$`+{\displaystyle \frac{2C_{\tau ,\tau }}{\rho _0}}{\displaystyle d^3p^{}\frac{f_\tau (\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}+{\displaystyle \frac{2C_{\tau ,\tau ^{}}}{\rho _0}}{\displaystyle d^3p^{}\frac{f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}`$
In the above $`\tau \tau ^{}`$ and $`\frac{\delta ^2}{\rho _n}=\frac{4\delta \rho _p}{\rho ^2}`$ and $`\frac{\delta ^2}{\rho _p}=\frac{4\delta \rho _n}{\rho ^2}`$.

The constants appearing in Eqs.(3.1) and (3.2) will be fixed by ensuring that properties of cold nuclear matter are reproduced. There $`f_\tau (\stackrel{}{r},\stackrel{}{p})=\frac{2}{h^3}\mathrm{\Theta }(p_f(\tau )p)`$. The integration in Eq.(3.1) is facilitated by noting that for a fixed $`\stackrel{}{p}(\stackrel{}{p}_1\stackrel{}{p}_2)/2`$, the centre of mass momentum can be integrated out to give

$`{\displaystyle _0^{p_f(\tau )}}{\displaystyle _0^{p_f(\tau ^{})}}d^3p_1d^3p_2g(\stackrel{}{p})`$ $`=`$ $`{\displaystyle _0^{q_f}}[{\displaystyle \frac{16\pi }{3}}(p_f^3(\tau )+p_f^3(\tau ^{}))8\pi p(p_f^2(\tau )+p_f^2(\tau ^{}))`$ (22)
$`+{\displaystyle \frac{16\pi }{3}}p^3{\displaystyle \frac{\pi }{p}}(p_f^2(\tau )p_f^2(\tau ^{}))^2]g(\stackrel{}{p})d^3p`$
where $`q_f=(p_f(\tau )+p_f(\tau ^{}))/2`$.

Fig.3 $`E_{sym}(\rho )`$ as a function of $`\rho `$ for the two choices of $`x`$; 0(left panel) and 1(right panel).

For completeness, for cold matter, we write down the values of the integrals appearing in Eqs.(3.1) and (3.2)

$`{\displaystyle d^3pd^3p^{}\frac{f_\tau (\stackrel{}{r},\stackrel{}{p})f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}`$ $`=`$ $`\left({\displaystyle \frac{2}{h^3}}\right)^2{\displaystyle \frac{4}{3}}\pi ^2\mathrm{\Lambda }^2[\{q_f{\displaystyle \frac{\mathrm{\Lambda }}{2}}\mathrm{arctan}({\displaystyle \frac{2q_f}{\mathrm{\Lambda }}})\}4(p_f^3(\tau )+p_f^3(\tau ^{}))`$ (25)
$`\{3(p_f^2(\tau )+p_f^2(\tau ^{}))+{\displaystyle \frac{\mathrm{\Lambda }^2}{2}}\}q_{f}^{}{}_{}{}^{2}+q_{f}^{}{}_{}{}^{4}`$
$`+\{{\displaystyle \frac{3\mathrm{\Lambda }^2}{4}}(p_f^2(\tau )+p_f^2(\tau ^{}))+{\displaystyle \frac{\mathrm{\Lambda }^4}{8}}{\displaystyle \frac{3}{8}}(p_f^2(\tau )p_f^2(\tau ^{}))^2\}ln(1+{\displaystyle \frac{4q_{f}^{}{}_{}{}^{2}}{\mathrm{\Lambda }^2}})]`$
Similarly, the value of the integral in Eq.(3.3) is

$`{\displaystyle d^3p^{}\frac{f_\tau (\stackrel{}{r},\stackrel{}{p}^{})}{1+(\stackrel{}{p}\stackrel{}{p}^{})^2/\mathrm{\Lambda }^2}}`$ $`=`$ $`{\displaystyle \frac{2}{h^3}}\pi \mathrm{\Lambda }^3[{\displaystyle \frac{p_f^2(\tau )+\mathrm{\Lambda }^2p^2}{2p\mathrm{\Lambda }}}ln{\displaystyle \frac{(p+p_f(\tau ))^2+\mathrm{\Lambda }^2}{(pp_f(\tau ))^2+\mathrm{\Lambda }^2}}`$ (27)
$`+{\displaystyle \frac{2p_f(\tau )}{\mathrm{\Lambda }}}2\{\mathrm{arctan}{\displaystyle \frac{p+p_f(\tau )}{\mathrm{\Lambda }}}\mathrm{arctan}{\displaystyle \frac{pp_f(\tau )}{\mathrm{\Lambda }}}\}]`$
Fig.4 One body potential in cold matter for $`\delta =0.4`$. As in fig. 3, here also the two panels are for $`x`$=0 and 1. Comparison with one-body potential as obtained with a Gogny interaction has also been shown. The solid line is the single particle potential for symmetric matter with MDI.

We fix the force parameters by first optimally reproducing the variation of $`U(\rho _0,\delta ,p)`$ with $`p`$ with that obtained by using a Gogny force (Brueckner-Hartee-Fock calculation gives similar results ). This fixes $`\mathrm{\Lambda }`$ of Eq.(3.3). Other parameters are then fixed by values of saturation density (0.16 $`fm^3`$), binding energy (16 MeV) at saturation density, compressibility ($``$ 210 MeV) of N=Z nuclear matter at saturation density and symmetry energy at $`N=Z`$ ($``$ 30 MeV). The values of the parameters for two choices (1) $`x`$=0 (denoted by MDI(0)) and (2) $`x`$=1 (denoted by MDI(1)) in Eq.(3.2) as in Gogny (this causes the value of the symmetry energy to bend over as a function of density) and are given in the table. Fig. 3 shows the behaviour of symmetry energy for the two choices. In fig. 4 we show the one body potential in cold matter for $`\delta =0.4`$.

Before we close this section we like to add that BUU calculations with momentum dependent interactions are only slightly more complicated in the asymmetric case compared to the case where no distinction is made between neutrons and protons. In numerical solutions of Boltzmann equations, phase space densities are simulated by test particles, characterised by a position and a momentum. Now there will also be a tag on their charge but no major addition to the codes  are needed.

## IV Reduction to the GBD form

Although we do recommend the full formalism of the above section be implemented, it is possible to reduce the above to a GBD (Gale, Bertsch and Das Gupta) form. A GBD potential, extended to asymmetric matter, already exists and is called BGBD . There the extensions were made such that in the $`N=Z`$ case one gets back exactly the original parameters . The parameter $`x`$ (eq.(3.1) was chosen to be 1/15. In contrast, here we have chosen the momentum dependence of the Gogny potential as a reference curve and chosen $`x`$ at 0 or 1. The values of the paramaters of GBD(0) ($`x`$=0), GBD(1) ($`x`$=1) and BGBD ($`x`$=1/15) are given in the table.

We write the potential energy density coming from the momentum dependent part as

$`V_{mom}(\rho ,\delta )={\displaystyle \frac{1}{\rho _0}}{\displaystyle \underset{\tau ,\tau ^{}}{}}C_{\tau ,\tau ^{}}\rho _\tau {\displaystyle \frac{f_\tau ^{}(\stackrel{}{r},\stackrel{}{p})d^3p}{1+(\stackrel{}{p}<\stackrel{}{p}>_\tau )^2/\mathrm{\Lambda }^2}}`$ (28)
The one-body potential generated by this piece of potential energy density for given $`\rho ,\delta ,\stackrel{}{p}`$ and a given $`\tau `$ is

$`U_{mom}(\rho ,\delta ,\stackrel{}{p},\tau )`$ $`=`$ $`{\displaystyle \frac{C_{\tau ,\tau }}{\rho _0}}\left[{\displaystyle \frac{f_\tau (\stackrel{}{r},\stackrel{}{p}^{})d^3p^{}}{1+(\stackrel{}{p}^{}<\stackrel{}{p}>_\tau )^2/\mathrm{\Lambda }^2}}+{\displaystyle \frac{\rho _\tau }{1+(\stackrel{}{p}<\stackrel{}{p}>_\tau )^2/\mathrm{\Lambda }^2}}\right]`$ (30)
$`+{\displaystyle \frac{C_{\tau ,\tau ^{}}}{\rho _0}}\left[{\displaystyle \frac{f_\tau ^{}(\stackrel{}{r},\stackrel{}{p}^{})d^3p^{}}{1+(\stackrel{}{p}^{}<\stackrel{}{p}>_\tau ^{})^2/\mathrm{\Lambda }^2}}+{\displaystyle \frac{\rho _\tau ^{}}{1+(\stackrel{}{p}<\stackrel{}{p}>_\tau ^{})^2/\mathrm{\Lambda }^2}}\right]`$
In the above equation, $`\tau ^{}\tau `$. As expected, the values of the constants in the force will have to be recalculated to reproduce the saturation properties. These are given in the table. The GBD potential $`U(\rho ,\delta ,p,\tau )`$ is plotted in Fig.5. This does not track the Gogny potential as faithfully as the more sophisticated version of section III does. The reader might wonder why the value of $`\mathrm{\Lambda }`$ is significantly bigger in GBD as opposed to in MDI(0) and MDI(1). We have tried to fit the variation of $`U`$ with $`k`$ as obtained in Gogny potential (or the Brueckner-Hartree-Fock calculation) by adjusting the value of $`\mathrm{\Lambda }`$. For $`p`$ close $`p_F`$, the contribution to $`U(p)`$ in section III comes mainly from $`p^{}`$ near $`p_F`$ whereas in GBD what counts is $`\stackrel{}{p}<\stackrel{}{p}>_{ave}`$ : since $`,\stackrel{}{<}p>_{ave}`$ is zero, one requires a different value of $`\mathrm{\Lambda }`$ to mimic the variation with $`p`$. This point was not appreciated in .

Fig.5 One body potential in cold matter for $`\delta =0.4`$. As in fig. 4, here also the two panels are for $`x`$=0 and 1. Comparison with one-body potential as obtained with a Gogny interaction has also been shown. The solid line represents the single particle potential for symmetric matter with the GBD potential.

## V Momentum dependence of the symmetry potential

The single particle potentials derived in the previous sections can be used directly in transport model calculations. They combine the density, momentum and isospin dependences of both the isoscalar and symmetry potentials in a nontrivial way. In this section we evaluate the strength of the momentum dependence of the symmetry potential. To the leading order in $`\delta `$, the single nucleon potential can be cast to the form

$$U_{n/p}(\rho ,\stackrel{}{p},\delta )U_0(\rho ,\stackrel{}{p})\pm U_{sym}(\rho ,\stackrel{}{p})\delta $$
(31)
in accordance with the Lane potential, where the $`\pm `$ sign is for neutrons and protons, respectively. Thus the symmetry potential can be evaluated from $`U_{sym}(\rho ,\stackrel{}{p})=(U_nU_p)/2\delta `$. Shown in Fig. 6 are the symmetry potentials as a function of $`k`$ for the three densities. It is seen that the $`U_{sym}(\rho ,\stackrel{}{p})`$ is strongly momentum dependent for $`k5fm^1`$ in all models considered. Moreover, this dependence is particularly stronger at high densities. By construction, the results for Gogny and MDI(1) are very close. By comparing the results with $`x=0`$ and $`x=1`$, it is seen that the momentum dependence is rather different mainly for $`\rho /\rho _0=2`$. This is because the symmetry energies with $`x=0`$ and $`x=1`$ are significantly different only in the region of $`\rho \rho _0`$ as shown in Fig. 3.

Fig.6 Dependence of $`(U_NU_P)/2\delta `$ on momentum for different potentials at different densities. All the figures are for $`\delta `$=0.4

To our best knowledge, the symmetry potential has been assumed to be momentum independent in all previous studies. Our results above indicate that this is only a good approximation at high momenta. At momenta less than about 1 GeV/c, the momentum dependence of the symmetry potential is important.

## VI Summary

The goal of this work was to obtain parameters of a momentum dependent mean field potential which is applicable to highly isospin asymmetric nuclear matter but easy enough to use in a transport model calculation for heavy ion collisions at intermediate to high energies. We used the Gogny interaction as a guide. Published results of Brueckner-Hartree-Fock calculations were also used for choosing parameters. We have two versions. They both are quite flexible in the sense that parameters can be easily chosen according to that experimental data on binding energy, compressibility etc. In each of these versions we have proposed two sets of parameters: one, in which the symmetry energy continues to rise as a function of density and another one where the symmetry energy first rises and then falls off with density. Our opinion is that introducing momentum dependence in symmetry potential will not not make transport model simulations much harder or longer than they already are. Implementations of these potentials in BUU transport models are in progress.

Transport models with the momentum dependent symmetry potentials are more reliable tools for investigating the density dependence of nuclear symmetry energy and thus the EOS of neutron-rich matter. Knowledge on the symmetry energy $`E_{sym}(\rho )`$ is essential for understanding not only the structure of radioactive nuclei but also many key issues in astrophysics. For instance, the $`E_{sym}(\rho )`$ determines uniquely the proton fraction in neutron stars at $`\beta `$ equilibrium. A continuously rising symmetry energy leads to a growing proton fraction with increasing density, thus allowing for the fast cooling of protoneutron stars through the direct URCA process. A falling symmetry energy at high densities forbides the direct URCA process to happen; moreover, it favors the formation of pure neutron domains in the cores of neutron stars. Nuclear reactions induced by neutron-rich nuclei provide a great opportunity to pin down the density dependence of nuclear symmetry energy. Among the experimental observables that have been found to be sensitive to the symmetry potential, the neutron/proton ratio of pre-equilibrium nucleon emissions, neutron-proton differential flow and correlation functions, as well as the proton elliptic flow at high transverse momenta are expected to be most sensitive to the momentum dependence of the symmetry potential. These observables will be studied with the improved BUU transport models using the momentum dependent symmetry potentials. These results will be reported in a forthcoming publication.

## VII Acknowledgement

S. Das Gupta acknowledges a very useful discussion with J.M. Pearson. Bao-An Li acknowledges the warm hospitality he received at McGill University where this work started. This research is supported in part by the Natural Sciences and Engineering Research Council of Canada and by the National Science Foundation of the United States under grant No. PHY-0088934.