text:
# Reversible coagulation-fragmentation processes and random combinatorial structures: asymptotics for the number of groups

## 1  Summary.

Our main result is a central limit theorem (Theorem $`4.6`$) for the number of groups at steady state for a class of reversible CFP’s and for the corresponding class of RCS’s.

In Section 2, we provide a definition of a reversible k-CFP admitting interactions of up to $`k`$ groups, as a generalization of the standard 2-CFP. The steady state of the processes considered is fully defined by a parameter function $`a0`$ on the set of integers. It was observed by F. Kelly (, p. 183) that for all $`2kN`$ the k-CFP’s have the same invariant measure on the set of partitions of a given integer $`N(=`$ the number of particles).

Section 3 explains the idea of A. Khintchine’s probabilistic method for derivation of asymptotic formulae. In the spirit of the method, we construct a representation of the probability function of the number of groups via the probability function of the sum of i.i.d. random variables.

In section 4 we study the case when the parameter function $`a`$ is of the form: $`a(k)=qk^{p1},k1,q,p>0.`$ We prove a local and a central normal limit theorems for the number of groups at equilibrium, as $`N\mathrm{}`$. To achieve this, we employ a new (for this field) tool: the Poisson summation formula. We conclude the section by providing some intuition for the main result of the paper.

In Section 5 we recall that the invariant measure of a reversible CFP can be viewed as a joint distribution of components of a RCS, known as an assembly. Comparing our main result with the known ones for RCS’s, we identify $`p=0`$ as a point of phase transition of the invariant measure, as $`N\mathrm{}`$. We also provide a few examples of RCS’s that conform to the case $`p>0`$.

## 2 CFP’s with multiple interactions: Definition.

Following , , we treat a CFP as a continuous-time Markov chain on the finite set $`\mathrm{\Omega }_N=\{\eta \}`$ of all partitions $`\eta =(n_1,\mathrm{},n_N)`$ of a given integer $`N`$:

$$\underset{j=1}{\overset{N}{}}jn_j=N,n_j0,j=1,\mathrm{},N.$$
(2.1)
Here $`N`$ codes the total population of indistinguishable particles partitioned into groups(=clusters) of different sizes. A group of size $`j2`$ may split into a number, say $`s,2sj`$, of groups of sizes $`j_1,\mathrm{},j_s:j_1+\mathrm{}+j_s=j`$ and, conversely, the above $`s`$ groups may coagulate into one large group of size $`j`$. We will call these $`s`$-interactions(=s-transitions), $`J_s`$ -fragmentation and $`J_s`$\- coagulation respectively, where $`J_s=(j_1,\mathrm{},j_s).`$ A stochastic process that admits interactions of up to $`kN`$ groups will be denoted k-CFP. Note that both types of the interactions conserve the total number of particles.

We now provide a formal definition of a k-CFP that naturally extends the definition of the standard 2-CFP (see , ). A k-CFP is given by the rates of infinitesimal transitions that are assumed to depend only on the sizes of interacting groups. For given $`N`$ and $`s,\mathrm{\hspace{0.33em}2}sk`$, let $`\psi _s,\varphi _s0`$ be a pair of functions defined on the same set

$$𝕁_s=\{J_s=(j_1,\mathrm{},j_s):j_l>0,l=1,\mathrm{},s,j_1+\mathrm{}j_sN\}$$

of s-tuples of positive integers, depicting the sizes of interacting groups. The functions $`\psi _s,\varphi _s`$ are the rates of $`J_s`$-coagulation and $`J_s`$-fragmentation respectively. Both functions are assumed to be invariant w.r.t. all $`s!`$ permutations of $`j_1,\mathrm{},j_s.`$ To complete the definition of a k-CFP, it is left to determine the total rates of all possible s-transitions from one partition $`\eta \mathrm{\Omega }_N`$ to another. Assume that the given $`J_s=(j_1,\mathrm{},j_s)𝕁_s`$ and $`\eta =(n_1,\mathrm{},n_N)\mathrm{\Omega }_N`$ are such that $`n_lm_l1,lJ_s`$, where $`m_l`$ counts the number of components in $`J_s`$ that are equal to $`l`$. In other words, we assume, that a $`J_s`$-coagulation of some groups in the partition $`\eta `$ is possible. Clearly, a given $`J_s`$-coagulation of any groups in $`\eta `$ transforms $`\eta `$ into the same partition that will be denoted $`\eta ^{(J_s)}\mathrm{\Omega }_N.`$ Similarly, the result of $`J_s`$-fragmentation of any groups in $`\eta \mathrm{\Omega }_N`$ will be denoted $`\eta _{(J_s)}\mathrm{\Omega }_N`$. By a simple combinatorial calculation, the total rate $`\mathrm{\Psi }_s(J_s;\eta )`$ of all possible $`J_s`$\- coagulations at the partition $`\eta `$ is

$$\mathrm{\Psi }_s(J_s;\eta )=\left(\underset{j_l}{}\frac{n_{j_l}!}{(n_{j_l}m_l)!}\right)\psi _s(J_s),J_s𝕁_s,$$
(2.2)
where the product is taken over all distinct components $`j_l`$ of $`J_s`$. By the same logic we define the total rate $`\mathrm{\Phi }(J_s;\eta )`$ of all possible $`J_s`$-fragmentations at the partition $`\eta `$ to be equal to

$$\mathrm{\Phi }_s(J_s;\eta )=n_{|J_s|}\varphi _s(J_s),J_s𝕁_s,$$
(2.3)
where $`|J_s|=j_1+\mathrm{}+j_s.`$ Now we see that a k-CFP is fully defined by the $`(k1)`$ pairs of functions $`\psi _s,\varphi _s,s=2,\mathrm{},k.`$

In this paper we will be concerned only with reversible CFP’s. Define the ratio of s-interactions

$$q_s(J_s)=\{\begin{array}{cc}\frac{\psi _s(J_s)}{\varphi _s(J_s)},\hfill & \text{if}J_s:\psi _s(J_s)\varphi _s(J_s)>0,J_s𝕁_s,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$
(2.4)
A natural extension of Theorem 1 in  gives the following characterization of reversible k-CFP’s with positive transition rates.

###### Proposition 2.1

Let

$$\psi _s(J_s)\varphi _s(J_s)>0,J_s𝕁_s,s=2,\mathrm{},k.$$
(2.5)
Then the corresponding k-CFP is reversible iff the ratios $`q_s,s=2,\mathrm{},k`$ are of the form:

$$q_s(J_s)=\frac{a(|J_s|)}{_{l=1}^sa(j_l)},J_s=(j_1,\mathrm{},j_s)𝕁_s,s=2,\mathrm{},k,$$
(2.6)
where $`a>0`$ is a function on the set of positive integers.

The proof is deferred until after Proposition 2.2.

A historical remark: It was already noted in , that the characterization of reversible 2-CFP’s was motivated by the following two completely independent lines of research: the seminal paper of F. Spitzer  on nearest-particle systems and F. Kelly’s and P. Whittle’s works in the 1970-s on networks and clustering process in polymerization (see , ). $`\mathrm{}`$

We will call $`a`$ the parameter function of a reversible CFP and we will write $`a_j=a(j),j1`$. The following result gives the explicit form of the steady state of the processes.

###### Proposition 2.2

(,,Ch.8)

For a given $`N`$, all k-CFP’s, $`k=2,\mathrm{},N`$ satisfying (LABEL:\[)03\] for some parameter function $`a`$ have the same invariant measure $`\mu _N`$ on $`\mathrm{\Omega }_N:`$

$$\mu _N(\eta )=(c_N)^1\frac{a_1^{n_1}a_2^{n_2}\mathrm{}a_N^{n_N}}{n_1!n_2!\mathrm{}n_N!},\eta =(n_1,\mathrm{},n_N)\mathrm{\Omega }_N,$$
(2.7)
where $`a_j>0,j=1,2,\mathrm{},N`$ and $`c_N=c_N(a_1,\mathrm{},a_N)`$ is the partition function of the measure $`\mu _N:`$

$$c_0=1,c_N=\underset{\eta \mathrm{\Omega }_N}{}\frac{a_1^{n_1}a_2^{n_2}\mathrm{}a_N^{n_N}}{n_1!n_2!\mathrm{}n_N!},N1.$$
(2.8)
Proof: We have to show that the measure $`\mu _N`$ given by $`(\text{LABEL:[})05]`$ satisfies the detailed balance condition:

$$\mu _N(\eta )V(\eta ,\xi )=\mu _N(\xi )V(\xi ,\eta ),\eta ,\xi \mathrm{\Omega }_N,$$
(2.9)
where $`V(\eta ,\xi )`$ is the total rate of the infinitesimal (in time) transition from $`\eta `$ to $`\xi `$. If $`V(\eta ,\xi )=0`$, then (LABEL:\[)401\] is trivially true. If now $`\xi =\eta ^{(J_s)}`$ for some $`J_s𝕁_s,\mathrm{\hspace{0.25em}2}sk`$, then we see from (LABEL:\[)05\] and (LABEL:\[)01\] that

$$\frac{\mu _N\left(\eta ^{(J_s)}\right)}{\mu _N(\eta )}=\frac{\mathrm{\Psi }_s(J_s;\eta )}{\mathrm{\Phi }_s(J_s;\eta ^{(J_s)})}.$$
(2.10)
In a similar manner one can verify (LABEL:\[)401\] in the case $`\xi =\eta _{(J_s)}`$ for some $`J_s𝕁_s`$. $`\mathrm{}`$
Now we are in position to give the
Proof of Proposition 2.1: If a k-CFP, $`k2`$ is reversible, then it follows from , Theorem 1, that (LABEL:\[)03\] should hold for s=2. The latter fact implies that the unique invariant measure of all k-CFP-s, $`k=2,\mathrm{},N`$ is given by (LABEL:\[)05\]. Consequently, (LABEL:\[)03\] should hold for $`3sk`$, by Proposition 2.2. The same reasoning proves the converse part of the claim. $`\mathrm{}`$
The preceding discussion shows that the steady state of a reversible k-CFP is uniquely determined by a parameter function $`a`$.

## 3 Khintchine’s type representation for the probability function of the number of groups.

Our objective will be the study of the asymptotic behavior, as $`N\mathrm{},`$ of the number of groups $`\nu _N`$ at equilibrium given by the measure $`\mu _N`$. It follows from (LABEL:\[)05\] and (LABEL:\[)08\] that

$$(\nu _N=n)=(c_N)^1\left(\underset{\eta \mathrm{\Omega }_N}{}\frac{a_1^{n_1}a_2^{n_2}\mathrm{}a_N^{n_N}}{n_1!n_2!\mathrm{}n_N!}\mathrm{𝟏}_{(_{k=1}^Nn_k=n)}\right),nN.$$
(3.11)
As in , , our tool will be the probabilistic method by A. Khintchine introduced in the 1950’s in his seminal book  . The idea of the method is to construct the representation of the quantity of interest via the probability function of a sum of independent integer valued random variables depending on a free parameter, and subsequent implementation of a local limit theorem. This allows for the derivation of the desired asymptotic formula.

In number theory, the implementation of Khintchine’s method was developed by G. Freiman *et al* (for references see ,). In particular, a general scheme of the method for asymptotic problems related to partitions was outlined by G. Freiman and J. Pitman in . The method was applied to CFP’s for the first time in , for derivation of the asymptotic formula for the partition function of the measure $`\mu _N`$ in the case $`a_kk^{p1},k\mathrm{},p>0.`$ In  the method was used for the study of the asymptotic behavior of some quantities related to clustering of groups at the steady state, when $`a_kk^{p1}L(k),k\mathrm{},p>0,`$ where $`L`$ is a slowly varying function.

Though the implementation of Khintchine’s method goes along the standard scheme, the related asymptotic analysis varies from problem to problem. In contrast to the aforementioned research, the problem treated in the present paper requires knowledge of the second term in the asymptotic expansions considered. In light of this, we employ the Poisson summation formula, a new tool for this field.

We will always assume that the parameter functions $`a`$ considered are positive and s.t. the power series

$$\underset{k=1}{\overset{\mathrm{}}{}}a_kx^k,x$$
(3.12)
has a finite radius of convergence $`R>0`$. Since the transformation $`a_kh^ka_k,h>0,k=1,2,\mathrm{},N`$ does not change the measure $`\mu _N`$ given by (LABEL:\[)05\], we assume w.l.g. that $`R=1.`$
Now let $`\xi _1,\mathrm{},\xi _n`$ be i.i.d. integer valued nonzero r.v.’s defined by

$$(\xi _1=l)=\frac{a_le^{\delta l}}{S(\delta )},a_l>0,l1,$$
(3.13)
where $`\delta >0`$ is a free parameter and

$$S(\delta )=\underset{k=1}{\overset{\mathrm{}}{}}a_ke^{\delta k}.$$
(3.14)
Note that the r.v. $`\xi _1`$ has finite moments of all orders for all $`\delta >0,`$ since $`R=1`$. We start with the following representation of the probability $`(\nu _N=n).`$

###### Lemma 3.3

Define

$$T_n=\underset{k=1}{\overset{n}{}}\xi _k,n1,$$

where $`\xi _k,k=1,2,\mathrm{}`$ are i.i.d. r.v.’s given by (LABEL:\[)1\],(LABEL:\[)2\]. Then

$`(\nu _N=n)=(c_Nn!)^1S^n(\delta )e^{\delta N}(T_n=N),\delta >0.`$ (3.15)
Proof: It follows from (LABEL:\[)1\],(LABEL:\[)2\] that

$$(T_n=N)=\frac{\left({\displaystyle \underset{\eta \mathrm{\Omega }_N}{}}{\displaystyle \frac{a_1^{n_1}a_2^{n_2}\mathrm{}a_N^{n_N}}{n_1!n_2!\mathrm{}n_N!}}\mathrm{𝟏}_{(_{k=1}^Nn_k=n)}\right)n!}{S^n(\delta )e^{\delta N}}.$$
(3.16)
By (LABEL:\[)05\] and (LABEL:\[)09\] this implies the claim (LABEL:\[)3\].

Remark. (LABEL:\[)3\] has a form of a typical representation in Khintchine’s method. It can be also viewed as a version of the representation formula for the total number of components in the generalized scheme of allocation (see , Lemma 1.3.3).

Our study will be heavily based on the assumption

$$\underset{k=1}{\overset{\mathrm{}}{}}ka_k=\mathrm{}.$$
(3.17)
Denote

$$M_1=M_1(n;\delta )=𝐄T_n=n𝐄\xi _1=nS^1(\delta )\underset{k=1}{\overset{\mathrm{}}{}}ka_ke^{\delta k},\delta >0$$
(3.18)
and choose the free parameter $`\delta `$ as a solution of the equation

$$M_1=N,nN.$$
(3.19)
Such a choice of the free parameter is typical for Khintchine’s method (, p.110) and is designed to make the probability $`(T_n=N)`$ in (LABEL:\[)3\] large, as $`n,N\mathrm{},nN.`$

###### Lemma 3.4

Under condition (LABEL:\[)0011\], the equation (LABEL:\[)5\] has a unique solution $`\delta =\delta _{n,N}`$ for all $`nN.`$

Proof: We first show that $`M_1`$ is decreasing in $`\delta >0:`$

$$M_1^{}(n;\delta )=n\frac{S(\delta )_{k=1}^{\mathrm{}}k^2a_ke^{\delta k}+\left(_{k=1}^{\mathrm{}}ka_ke^{\delta k}\right)^2}{S^2(\delta )}<0,\delta >0,$$
(3.20)
by the Cauchy-Schwarz inequality. Consequently, $`M_1(n;0)=\underset{\delta >0}{sup}M_1(n;\delta ):=nM_1^{}`$ where $`M_1^{}`$ does not depend on $`n`$. If now $`S(0)<\mathrm{}`$, then (LABEL:\[)0011\] immediately implies $`M_1^{}=\mathrm{}.`$ In the case $`S(0)=\mathrm{}`$, supposing $`M_1^{}<\mathrm{}`$ leads to the contradiction:

$`0`$ $`>`$ $`M_1(n;\delta )nM_1^{}=n{\displaystyle \frac{_{k=1}^{\mathrm{}}a_ke^{\delta k}(kM_1^{})}{S(\delta )}}`$ (3.21)
$``$ $`nS^1(\delta )\left(A(\delta )+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k2M_1^{}}{}}ka_ke^{\delta k}\right),\delta >0,`$
where $`A(\delta )`$ is bounded for any $`\delta 0`$, while the sum in the brackets tends to $`+\mathrm{}`$ as $`\delta 0^+`$, by (LABEL:\[)0011\]. Hence, $`M_1(n;0)=+\mathrm{}`$. Finally,

$$M_1(n;\mathrm{})=n\underset{\delta \mathrm{}}{lim}\frac{_{k=1}^{\mathrm{}}ka_ke^{\delta (k1)}}{_{k=1}^{\mathrm{}}a_ke^{\delta (k1)}}=n.$$
(3.22)
## 4 A central limit theorem for the number of groups.

Our paper is devoted exclusively to the case when the parameter function $`a`$ has a polynomial rate of growth, namely:

$$a_k=qk^{p1},q,p>0,k1.$$
(4.23)
We first consider the case $`q=1.`$ The following lemma which is basic for our subsequent asymptotic analysis is a particular case of the Poisson summation formula (see ,).

###### Lemma 4.5

(, p.82). Let $`p>1`$ and $`Re(z)>0.`$ Then we have

$$\underset{k=1}{\overset{\mathrm{}}{}}e^{zk}k^{p1}=\mathrm{\Gamma }(p)\left(\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}(z+2\pi il)^p\right).$$
(4.24)
With the help of this remarkable identity we derive the following asymptotic formula that holds for $`p>0`$.

###### Lemma 4.6

If $`p>0,|z|0,Re(z)0^+,`$ then

$$\underset{k=1}{\overset{\mathrm{}}{}}e^{zk}k^{p1}=\mathrm{\Gamma }(p)\left(z^p+A(p)\right)+O(z),$$
(4.25)
where $`A(p)`$ is a constant which in the case $`p>1`$ is given explicitly by

$$A(p)=2(2\pi )^p\zeta (p)\mathrm{cos}\frac{\pi p}{2}.$$
(4.26)
(Here $`\zeta (p)`$ is the Riemann zeta function).

Proof: First consider the case $`p>1.`$ Let

$$G(z)=\underset{l=\mathrm{}}{\overset{1}{}}(z+2\pi il)^p+\underset{l=1}{\overset{\mathrm{}}{}}(z+2\pi il)^p,Re(z)0,p>1.$$
(4.27)
The two series in the RHS of (LABEL:\[)12a\] converge absolutely when $`Re(z)>0,`$ while, by straightforward calculations, $`G(0)=A(p),G^{}(0)=pA(p+1)`$ with $`A(p)`$ given by (LABEL:\[)12\]. Consequently, by (LABEL:\[)10\] we have for $`p>1`$

$`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}e^{zk}k^{p1}\mathrm{\Gamma }(p)z^p`$ $`=`$ $`\mathrm{\Gamma }(p)G(z)=\mathrm{\Gamma }(p)\left(G(0)+G^{}(0)z+o(z)\right),`$ (4.28)
$`|z|0,Re(z)0^+.`$
This proves (LABEL:\[)11\] for $`p>1.`$ So, we write for $`p>0`$

$$\underset{k=1}{\overset{\mathrm{}}{}}e^{zk}k^p=\mathrm{\Gamma }(p+1)\left(z^{p1}+A(p+1)\right)+O(z),|z|0,Re(z)0^+.$$
(4.29)
Next, integrating (LABEL:\[)12c\] w.r.t. $`z`$ gives (LABEL:\[)11\] with a constant $`A(p)`$ that is not known explicitly.

Now we are in a position to derive the asymptotic formula for the free parameter $`\delta `$.

###### Proposition 4.7

Assume that $`n`$ is s.t. $`\alpha :=\frac{N}{n}\mathrm{},`$ as $`N\mathrm{}`$. Then

$$\delta =p\alpha ^1\left(1A(p)p^p\alpha ^p\right)+o\left(\alpha ^{p1}\right),\alpha \mathrm{}.$$
(4.30)
Proof: First observe that in the case of the function $`a`$ considered, it follows from (LABEL:\[)4a\],(LABEL:\[)5\] and the proof of Lemma 3.2, that $`\alpha =\frac{N}{n}\mathrm{}`$ implies $`\delta _{n,N}0^+`$. Implementing (LABEL:\[)11\] gives

$`𝐄\xi _1`$ $`=`$ $`{\displaystyle \frac{{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k^pe^{\delta k}}{{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}k^{p1}e^{\delta k}}}=p{\displaystyle \frac{\delta ^{p1}+A(p+1)+O(\delta )}{\delta ^p+A(p)+O(\delta )}}=`$ (4.31)
$`=`$ $`p\delta ^1\left(1A(p)\delta ^p+o(\delta ^p)\right),\delta 0^+,p>0.`$
Now (LABEL:\[)5\] leads to

$$\delta =p\alpha ^1\left(1A(p)\delta ^p+o(\delta ^p)\right),\delta =\delta _{n,N}0^+,\alpha \mathrm{}.$$
(4.32)
Iterating this equation w.r.t. $`\delta `$ gives (LABEL:\[)13\].

We will focus now on asymptotics, as $`\alpha \mathrm{},`$ of the probability $`(T_n=N)`$ under $`\delta `$ given by (LABEL:\[)13\]. We have

$$(T_n=N)=(2\pi )^1\left(_\pi ^\pi \phi (t)e^{itN}𝑑t\right),$$
(4.33)
where $`\phi `$ is the characteristic function of $`T_n.`$

The basic idea of the Khintchine’s method is that for a wide class of models, choosing the free parameter from the condition (LABEL:\[)5\] guarantees that the main contribution to the integral in the RHS of (LABEL:\[)14\] comes from a set which is some neighborhood of zero.

We will demonstrate that this is in effect true in the case considered and will prove the normal local limit theorem for the sum $`T_n,`$ as $`\alpha \mathrm{}.`$ As a preliminary step, we verify the validity of the Lyapunov’s condition

$$\underset{n\mathrm{}}{lim}\frac{M_3}{M_2^{\frac{3}{2}}}=0,$$
(4.34)
where $`M_2=M_2(n;N),M_3=M_3(n;N)`$ are correspondingly the variance and the third central moment of $`T_n`$ under $`\delta `$ given by (LABEL:\[)13\].

We have

$$M_2=VarT_n=nVar\xi _1=n\left(S^1(\delta )\underset{k=1}{\overset{\mathrm{}}{}}k^2a_ke^{\delta k}\alpha ^2\right),\delta >0$$
(4.35)
and

$$M_3=E(T_nN)^3=nE(\xi _1\alpha )^3=n(E\xi _1^33\alpha E\xi _1^2+2\alpha ^3).$$
(4.36)
Applying now (LABEL:\[)11\] and (LABEL:\[)13\] gives the following asymptotic expressions for the moments considered:

$$M_2n\left(p(p+1)\delta ^2p^2\delta ^2\right)=np\delta ^2,$$

$$M_3n\left(p(p+1)(p+2)\delta ^33p^2(p+1)\delta ^3+2p^3\delta ^3\right)=2p\delta ^3n,n,\alpha \mathrm{}.$$
(4.37)
This proves (LABEL:\[)4f\]. The condition (LABEL:\[)4f\] provides the existence of $`\beta =\beta (n,\delta )>0`$ s.t.

$$\beta ^2M_2\mathrm{},\beta ^3M_30,n,\alpha \mathrm{}.$$
(4.38)
Explicitly, in view of (LABEL:\[)4c\],

$$\beta =\delta n^{\frac{1ϵ}{2}},0<ϵ<1/3$$
(4.39)
satisfies (LABEL:\[)4g\]. As it will be shown below, $`[\beta ,\beta ]`$ is just the required neighborhood of zero.

###### Lemma 4.8

(The local limit theorem for $`T_n`$.)

Set

$$n=Q_pN^{\frac{p}{p+1}}+sN^{\frac{p}{2p+2}},Q_p=p^1(\mathrm{\Gamma }(p+1))^{\frac{1}{p+1}},sR,p>0.$$
(4.40)
Then

$$(T_n=N)\frac{1}{\sqrt{2\pi VarT_n}},N\mathrm{}.$$
(4.41)
Proof: We write

$$I=I_1+I_2,$$
(4.42)
where $`I={\displaystyle _\pi ^\pi }\phi (t)e^{itN}𝑑t`$ and $`I_1,I_2`$ are integrals of the integrand of $`I`$ taken over the sets $`[\beta ,\beta ]`$ and $`[\pi ,\beta ][\beta ,\pi ]`$ respectively.

Step 1. We find the asymptotics of the integral $`I_1,`$ when $`\beta `$ is as given by (LABEL:\[)4n\].

By the definition of $`\alpha ,`$

$$\phi (t)e^{itN}=\phi _1^n(t),tR,$$
(4.43)
where $`\phi _1`$ is the characteristic function of the r.v. $`\xi _1\alpha .`$ Next, denoting by $`\phi _2`$ the characteristic function of the r.v. $`\left(p^{1/2}\delta \right)(\xi _1\alpha ),`$ we get from (LABEL:\[)4\]–(LABEL:\[)4c\]

$$\phi _2(t)=1\frac{1}{2}t^2+O(t^3),t0.$$
(4.44)
Combining this with the relationship

$$\phi _2(t)=\phi _1(p^{1/2}\delta t),tR,$$
(4.45)
(LABEL:\[)18\] becomes

$$\phi (t)e^{itN}=\left(1\frac{1}{2}(\sqrt{p}\delta ^1t)^2+O(\delta ^3t^3)\right)^n$$

$$\mathrm{exp}\left(\frac{1}{2}(\sqrt{np}\delta ^1t)^2+nO(\delta ^3t^3)\right),t\delta ^10,n\mathrm{}.$$
(4.46)
Consequently, by virtue of (LABEL:\[)4c\] and (LABEL:\[)4g\],

$$I_1\sqrt{\frac{2\pi }{VarT_n}},n,\alpha \mathrm{},p>0.$$
(4.47)
Step 2. We are to show that

$$I_2=o(I_1),N\mathrm{},p>0.$$
(4.48)
We apply (LABEL:\[)10\] with $`z=\delta it,t[\beta ,\pi ]`$ to obtain from (LABEL:\[)12a\] the following analog of (LABEL:\[)12b\]:

$$\underset{k=1}{\overset{\mathrm{}}{}}e^{zk}k^{p1}=\mathrm{\Gamma }(p)\left(z^p+G(it)+\delta G^{}(it)\right)+o(\delta ),p>1,t[\beta ,\pi ],\delta 0^+,$$
(4.49)
where $`|G(it)|,|G^{}(it)|const:=B(p),t[\beta ,\pi ],p>1.`$ By the same argument as for (LABEL:\[)12c\] the latter yields

$$|\phi _1(t)|\frac{(\delta ^2+\beta ^2)^{\frac{p}{2}}+B(p)+O(\delta )}{\delta ^p+A(p)+O(\delta )},t[\beta ,\pi ],p>0,\delta 0^+.$$
(4.50)
In the rest of this section it is always assumed that $`n`$ is specified as in (LABEL:\[)40\]. In view of (LABEL:\[)40\] we have

$$n\alpha ^p=\frac{n^{p+1}}{N^p}=p^p\mathrm{\Gamma }(p)+ϵ_N,p>0,N\mathrm{},$$
(4.51)
where

$$ϵ_N=(p+1)Q_p^psN^{\frac{p}{2p+2}}+\frac{p(p+1)}{2}Q_p^{p1}s^2N^{\frac{p}{p+1}}+O(N^{\frac{3p}{2p+2}}),p>0,N\mathrm{}.$$
(4.52)
This fact will be repeatedly used in our subsequent asymptotic analysis. Consequently,

$$n\delta ^p=\mathrm{\Gamma }(p)+p^pϵ_N+O\left(\alpha ^p\right),p>0,N\mathrm{}.$$
(4.53)
Employing (LABEL:\[)4200\] we further obtain from (LABEL:\[)223\] and (LABEL:\[)4n\]

$$|\phi _1(t)|^nO\left(\mathrm{exp}(\frac{p}{2}n^ϵ)\right),t[\beta ,\pi ],p>0,n\mathrm{},$$
(4.54)
for any $`0<ϵ<\frac{1}{3}`$. This together with (LABEL:\[)22\] proves (LABEL:\[)0101\].

To establish our main result, it is left to find the asymptotic formulae for the rest of the factors in (LABEL:\[)3\]. First we make use of the following result of :

$$c_N\frac{1}{\sqrt{2\pi B_N^2}}\mathrm{exp}\left(N\sigma +\underset{k=1}{\overset{N}{}}k^{p1}e^{k\sigma }\right),N\mathrm{},p>0,$$
(4.55)
where

$$\sigma =\sigma _N\left(\frac{N}{\mathrm{\Gamma }(p+1)}\right)^{\frac{1}{p+1}},p>0,N\mathrm{}$$
(4.56)
is the unique solution of the equation

$$\underset{k=1}{\overset{N}{}}e^{k\sigma }k^p=N,p>0,$$
(4.57)
while

$$B_N^2=\underset{k=1}{\overset{N}{}}e^{k\sigma }k^{p+1},p>0.$$
(4.58)
For our purpose, we need to know the second term in the asymptotic expansion (LABEL:\[)27\] of $`\sigma `$. In what follows we denote by $`ϵ_N`$ different quantities tending to zero, as $`N\mathrm{}.`$ By our asymptotic formula (LABEL:\[)11\],

$$\underset{k=1}{\overset{\mathrm{}}{}}k^pe^{k\sigma }=\mathrm{\Gamma }(p+1)\left(\sigma ^{(p+1)}+A(p+1)\right)+O(\sigma ),p>0,N\mathrm{},$$
(4.59)
where $`A(p)`$ is as in Lemma 4.2, while, by the Euler summation formula, we have for any $`\gamma >0,`$ and $`\sigma =\sigma _N`$ given by (LABEL:\[)27\],

$$\underset{k=N+1}{\overset{\mathrm{}}{}}k^pe^{k\sigma }_{N+1}^{\mathrm{}}x^pe^{\sigma x}𝑑x=\sigma ^{(p+1)}_{\sigma (N+1)}^{\mathrm{}}x^pe^x𝑑x=o(N^\gamma ),p>0,N\mathrm{}.$$
(4.60)
In view of (LABEL:\[)28\],(LABEL:\[)29\], the equation (LABEL:\[)25\] can be rewritten now as

$$\mathrm{\Gamma }(p+1)\left(\sigma ^{(p+1)}+A(p+1)\right)+O(\sigma )=N.$$
(4.61)
Consequently, we get

$$\sigma =\left(\frac{N}{\mathrm{\Gamma }(p+1)}A(p+1)+ϵ_N\right)^{\frac{1}{p+1}}=$$

$$\left(\mathrm{\Gamma }(p+1)\right)^{\frac{1}{p+1}}N^{\frac{1}{p+1}}+\frac{A(p+1)}{p+1}\left(\frac{\mathrm{\Gamma }(p+1)}{N}\right)^{\frac{p+2}{p+1}}+o(N^{\frac{p+2}{p+1}}),p>0,N\mathrm{}.$$
(4.62)
This yields

$$N\sigma =N^{\frac{p}{p+1}}(\mathrm{\Gamma }(p+1))^{\frac{1}{p+1}}+ϵ_N,p>0,N\mathrm{}.$$
(4.63)
Next, we obtain from (LABEL:\[)11\],(LABEL:\[)26\],(LABEL:\[)29\] and (LABEL:\[)30\]

$$B_N^2\left(\frac{N}{\mathrm{\Gamma }(p+1)}\right)^{\frac{p+2}{p+1}}\mathrm{\Gamma }(p+2),p>0,N\mathrm{}$$
(4.64)
and

$$\underset{k=1}{\overset{N}{}}k^{p1}e^{k\sigma }=\underset{k=1}{\overset{\mathrm{}}{}}k^{p1}e^{k\sigma }+ϵ_N=$$

$$\mathrm{\Gamma }(p)\left(\frac{N}{\mathrm{\Gamma }(p+1)}\right)^{\frac{p}{p+1}}+\mathrm{\Gamma }(p)A(p)+ϵ_N,p>0,N\mathrm{}.$$
(4.65)
Substituting the above expressions in (LABEL:\[)2444\] gives the desired asymptotic formula for $`c_N`$:

$$c_Nh_1N^{\frac{p+2}{2(p+1)}}\mathrm{exp}(h_2N^{\frac{p}{p+1}}+h_3),p>0,N\mathrm{},$$
(4.66)
where the constants $`h_i,i=1,2,3`$ are given by

$$h_1=\frac{\left(\mathrm{\Gamma }(p+1)\right)^{\frac{1}{2(p+1)}}}{\sqrt{2\pi (p+1)}},$$

$$h_2=(p+1)Q_p,$$

$$h_3=\mathrm{\Gamma }(p)A(p),p>0.$$
(4.67)
Next, (LABEL:\[)4111\], (LABEL:\[)40\] and (LABEL:\[)13\] give

$$\mathrm{exp}(\delta N)\mathrm{exp}\left(pnA(p)\mathrm{\Gamma }(p+1)\right),p>0,N\mathrm{},$$
(4.68)
while (LABEL:\[)4200\] gives

$`S^n(\delta )`$ $`=`$ $`\left(\mathrm{\Gamma }(p)\delta ^p+A(p)\mathrm{\Gamma }(p)+O(\delta )\right)^n\left(\mathrm{\Gamma }(p)\right)^n\delta ^{pn}\mathrm{exp}\left(A(p)\mathrm{\Gamma }(p)\right),`$ (4.69)
$`p>0,N\mathrm{}.`$
We again apply (LABEL:\[)4111\] to get

$$\delta ^{pn}p^{pn}\alpha ^{pn}\mathrm{exp}\left(A(p)\mathrm{\Gamma }(p+1)\right),$$

$$\alpha ^{pn}\frac{n^np^{pn}\mathrm{exp}\left(n\left(\frac{p^p}{\mathrm{\Gamma }(p)}ϵ_N\frac{1}{2}\left(\frac{p^p}{\mathrm{\Gamma }(p)}\right)^2ϵ_N^2+O(ϵ_N^3)\right)\right)}{\left(\mathrm{\Gamma }(p)\right)^n},p>0,N\mathrm{},$$
(4.70)
where $`ϵ_N`$ is given by (LABEL:\[)400\]. Observing that $`nϵ_N^30,N\mathrm{},`$ we write out now the asymptotic expressions for $`nϵ_N`$ and $`nϵ_N^2`$ to obtain

$$S^n(\delta )n^n\mathrm{exp}\left(\frac{s^2}{2d_p}s(p+1)N^{\frac{p}{2p+2}}+A(p)\left(\mathrm{\Gamma }(p)+\mathrm{\Gamma }(p+1)\right)\right),$$

$$d_p=\frac{Q_p}{p+1},p>0,N\mathrm{}.$$
(4.71)
Finally, substituting in (LABEL:\[)3\] the preceding asymptotic expansions and employing Stirling’s asymptotic formula gives the following result.

###### Theorem 4.9

(The local limit theorem for $`\nu _N`$).

Let $`n`$ be given as in (LABEL:\[)40\]. Then

$$(\nu _N=n)\frac{1}{\sqrt{2\pi d_p}}N^{\frac{p}{2p+2}}\mathrm{exp}\left(\frac{s^2}{2d_p}\right):=f(N;s),sR,p>0,N\mathrm{}.$$
(4.72)
This leads to our main result that says that almost all the mass of the probability distribution of the r.v. $`\nu _N`$ is concentrated , as $`N\mathrm{},`$ in a neighborhood of size $`O(N^{\frac{p}{2p+2}})`$ of the point $`Q_pN^{\frac{p}{p+1}}`$.

###### Theorem 4.10

(The central limit theorem for $`\nu _N`$).

$$\frac{\nu _NQ_pN^{\frac{p}{p+1}}}{\sqrt{d_p}N^{\frac{p}{2p+2}}}N(0,1),p>0,N\mathrm{},$$
(4.73)
where $``$ denotes weak convergence and $`d_p,Q_p`$ are as in (LABEL:\[)46\],(LABEL:\[)40\] respectively.

Proof: We provide a sketch of the proof that is done by the implementation of the standard technique of passing from the local theorem to the integral theorem (for more details see (, p.59, , p.81).

It follows from (LABEL:\[)47\] that for any $`ab,a,bR,`$

$$\left(\frac{\nu _NQ_pN^{\frac{p}{p+1}}}{N^{\frac{p}{2p+2}}}[a,b]\right)=\underset{sR_N}{}f(N;s),p>0,N\mathrm{},$$
(4.74)
where $`f(N;s)`$ is given by (LABEL:\[)47\] and

$$R_N=\{s[a,b]:n=Q_pN^{\frac{p}{p+1}}+sN^{\frac{p}{2p+2}}\}.$$

Since $`|R_N|=O(N^{\frac{p}{2p+2}}),`$ as $`N\mathrm{},`$ we have

$$\underset{sR_N}{}exp(\frac{s^2}{2d_p})N^{\frac{p}{2p+2}}_a^b\mathrm{exp}(\frac{x^2}{2d_p})dx,p>0,N\mathrm{}.\text{ }$$
(4.75)
Remark : Since the transformation $`a_kh^ka_k,h>0,k1`$ of the parameter function $`a`$ does not change the measure $`\mu _N`$, the results of our paper are true for $`a_k=h^kk^{p1},k1,h>0,p>0`$. $`\mathrm{}`$

One more extension of Theorem 4.6 is provided by the following result.

###### Theorem 4.11

: If $`\stackrel{~}{a_k}=qk^{p1},p>0,k1`$, where $`q>0`$ is a constant, then

$`{\displaystyle \frac{\nu _N\stackrel{~}{Q_p}N^{\frac{p}{p+1}}}{\sqrt{\stackrel{~}{d_p}}N^{\frac{p}{2p+2}}}}N(0,1),p>0,N\mathrm{},`$ (4.76)
where $`\stackrel{~}{Q_p}=q^{\frac{1}{p+1}}Q_p`$, $`\stackrel{~}{d_p}=q^{\frac{1}{p+1}}d_p`$.

Proof: Denote by $`\stackrel{~}{}`$ the quantities related to the parameter function $`\stackrel{~}{a}`$. We see from (LABEL:\[)1\],(LABEL:\[)25\] and (LABEL:\[)30\] that the r.v. $`\stackrel{~}{\xi _1}`$ has the same distribution as $`\xi _1`$ and that

$$\stackrel{~}{\sigma _N}q^{\frac{1}{p+1}}\sigma _N,N\mathrm{}.$$
(4.77)
Repeating the preceding asymptotic analysis gives the claimed change in the scaling induced by $`q>0`$. $`\mathrm{}`$
We conclude this section by providing some intuition for the scaling in the central limit Theorem 4.6. For this purpose we employ the following result established in  for the model in question. Denote by $`\overline{q}(\eta ),\underset{¯}{q}(\eta )`$ the size of the largest (resp. smallest) group in a random partition $`\eta \mathrm{\Omega }_N.`$ Then

$`\underset{N\mathrm{}}{lim}\left(N^{\frac{1}{p+1}ϵ}<\overline{q}(\eta )<N^{\frac{1}{p+1}+ϵ}\right)=1,p>0,`$ (4.78)
for all $`ϵ>0`$, while

$$\underset{N\mathrm{}}{lim}\left(\underset{¯}{q}(\eta )l\right)=\{\begin{array}{cc}0,\hfill & \text{if}l=N^\beta ,0<\beta 1,\hfill \\ \mathrm{exp}\left(\underset{j=1}{\overset{l1}{}}a_j\right),\hfill & \text{if}l2\text{is a fixed number}.\hfill \end{array}$$
(4.79)
From (LABEL:\[)60a\] and (LABEL:\[)61b\] one may conjecture, that for large $`N`$, the main ”mass”, of size $`O(N)`$, is partitioned into groups(=clusters) of sizes $`O\left(N^{\frac{1}{p+1}}\right)`$, while the rest of the mass, of size $`o(N)`$, is partitioned into groups of small sizes. Adopting this conjecture gives the expectation of the number of groups $`\nu _N`$ as $`O\left(\frac{N}{N^{\frac{1}{p+1}}}\right)=O\left(N^{\frac{p}{p+1}}\right)`$. Concerning the asymptotic variance of $`\nu _N`$, we note that the relationship $`Var(\nu _N)=O(𝔼\nu _N),N\mathrm{}`$ holds also for permutations, the Ewens sampling formula and for some other RCS’s (see , ).

## 5 RCS’s: examples and comparison with known results.

It was already observed in  and  that particular cases of the expression (LABEL:\[)05\] for the equilibrium measure $`\mu _N`$ conform to joint distributions of components of a variety of RCS’s, known as assemblies.

By a combinatorial structure (CS) of a size $`N`$ one means a union of nondecomposable components(=components) of different sizes. Formally, such a structure is given by the two sets of integers $`\{p_N,N1\}`$ and $`\{m_N,N1\}`$ that count respectively the total number of instances of size $`N`$ and the number of components of size $`N.`$ An example of a CS is a graph on $`N`$ vertices treated as a union of its connected components. Therefore, an instance of a CS of size $`N`$ is given by a partition $`\eta =(n_1,\mathrm{},n_N)\mathrm{\Omega }_N`$, where $`n_k`$ is the number of components of size $`k`$ in the instance. By assuming that for a given $`N`$ an instance is chosen randomly from all $`p_N`$ instances, one induces a RCS that is completely determined by the random component size counting process, the latter being a random vector with values in $`\mathrm{\Omega }_N.`$ With an obvious abuse of notation we denote the random vector by $`\eta .`$ A remarkable fact in the theory of RCS’s is that a great variety of them obey the conditioning relation

$$(\eta )=(Z_1,\mathrm{},Z_N|\underset{j=1}{\overset{N}{}}jZ_j=N),$$
(5.1)
where $`Z_1,\mathrm{},Z_N`$ are independent integer valued r.v.’s. In particular, a class of CS’s known as assemblies is characterized by the fact that $`Z_k`$ are Poisson r.v.’s s.t. $`Z_kPo(m_k/k!),k=1,2,\mathrm{},N.`$ Consequently, by a straightforward calculation we find from (LABEL:\[)90\] that for assemblies

$$(\eta )=\mu _N(\eta ),\eta \mathrm{\Omega }_N,$$
(5.2)
with $`a_k=\frac{m_k}{k!}`$ and $`c_k=\frac{p_k}{k!},k=1,2,\mathrm{},N`$. In particular, the case $`m_k\theta (k1)!,\theta >0`$ conforms to logarithmic RCS’s that encompass permutations ($`m_k=(k1)!`$) and the Ewens sampling formula ($`m_k=\theta (k1)!,\theta >0)`$. The novel theory of general logarithmic RCS’s is presented in . From the analytical point of view, the common feature of logarithmic RCS’s is that they do not obey the condition $`p>0`$ adopted in the present paper. As a result, our asymptotic analysis based on the Poisson summation formula is not applicable for logarithmic RCS’s.

The case of permutations has a long history. For this case, the seminal result by V.L. Goncharov (1944) and L. Shepp and S. Lloyd (1966) states the following central limit theorem for $`\nu _N`$:

$$\frac{\nu _N\mathrm{log}N}{\sqrt{\mathrm{log}N}}N(0,1).$$
(5.3)
A version of (LABEL:\[)50\] for general logarithmic RCS’s is also known.
L. Mutafchiev () proved a local limit theorem for $`\nu _N`$ under some assumptions on the asymptotic behavior as $`x1`$ of the generating function of the sequence $`\{c_n\}_1^{\mathrm{}}`$ . It can be verified that in the case $`p>0`$ the assumptions $`(2.7)`$ and $`(2.8)`$ in  do not hold. However, Mutafchiev (,p.425) conjectured, that a result similar to his Theorem 2.4 holds for a wider class of RCS’s. Our Theorem 4.6 confirms this conjecture. In fact, in the case $`p>0`$, we have $`𝔼\nu _NQ_pN^{\frac{p}{p+1}}`$ and it is not hard to see that, in agreement with the claim in , $`Q_pN^{\frac{p}{p+1}}S(e^{\sigma _N}),N\mathrm{}`$, where $`\sigma _N`$ is as in (LABEL:\[)27\].
A. Barbour and B. Granovsky () explored the case when $`Z_k`$ in (LABEL:\[)90\] are quite general r.v.’s obeying $`𝔼Z_kk^{p1},k\mathrm{},p<0.`$ (This includes assemblies with $`m_kk^{p1}k!,k\mathrm{},p<0`$). It was shown in () that such RCS’s exhibit a completely different asymptotic behavior. In particular, in this case $`\nu _N`$ is finite with probability $`1`$.
Hence, comparing the asymptotic behavior of $`\nu _N`$ (e.g. $`𝔼\nu _N)`$ in the cases $`p<0`$ and $`p0`$, one sees that $`p=0`$ can be viewed as a point of phase transition for the measure $`\mu _N`$ as $`N\mathrm{}`$.

In conclusion, we provide a few examples of assemblies that conform to the setting of the present paper: Example 1. Forests of labelled and colored linear trees. A linear tree (see ) is a graph with no cycles, where each vertex has no more than two neighbors. Assuming that a vertex is labelled and is colored into one of $`q(q1)`$ colors gives $`m_k=q^kk!`$, $`a_k=q^k,k1`$. So, by the remark following Theorem 4.6, this RCS corresponds to the case $`p=1`$.

Example 2. Forests of labelled rooted linear trees. In this case we have $`m_k=kk!,k1`$, which gives $`a_k=k,k1`$, that corresponds to the case $`p=2`$.

Example 3. Compositions. (see ) Consider a space $`\mathrm{{\rm Y}}_N`$ of ordered m-tuples$`𝐱=(x_1,\mathrm{},x_m),m=1,\mathrm{},N`$, where $`x_i`$ are positive integers, summing to $`N`$. In other words, $`\mathrm{{\rm Y}}_N`$ is a space of all ordered partitions (=compositions) of $`N`$. We define the probability measure $`\lambda _N`$ on $`\mathrm{{\rm Y}}_N`$:

$$\lambda _N(𝐱)=\frac{1}{(m(𝐱))!}(r_N)^1,𝐱=(x_1,\mathrm{},x_m)\mathrm{{\rm Y}}_N,$$
(5.4)
where $`(r_N)^1`$ is the normalizing constant and $`m=m(𝐱)`$ is the number of components of $`𝐱`$. This means that $`\lambda _N`$ prescribes to $`𝐱𝚼_𝐍`$ the weight $`\frac{1}{(m(𝐱))!}`$. Denote $`n_i=n_i(𝐱),i=1,\mathrm{},N`$ the number of components of $`𝐱𝚼_𝐍`$ that are equal to $`i`$. Then we have

$`r_N={\displaystyle \underset{m=1}{\overset{N}{}}}{\displaystyle \frac{1}{m!}}{\displaystyle \underset{𝐱𝚼_𝐍:𝐦(𝐱)=𝐦}{}}1={\displaystyle \underset{m=1}{\overset{N}{}}}{\displaystyle \frac{1}{m!}}{\displaystyle \underset{\eta \mathrm{\Omega }_N:|\eta |=m}{}}{\displaystyle \frac{m!}{n_1!\mathrm{}n_m!}}={\displaystyle \underset{\eta \mathrm{\Omega }_N}{}}{\displaystyle \frac{1}{n_1!\mathrm{}n_N!}},`$ (5.5)
where, as before, $`\mathrm{\Omega }_N`$ is the set of all (unordered) partitions of $`N`$, while $`|\eta |=n_1+\mathrm{}+n_N`$. (LABEL:\[)999\] says that $`r_N=c_N`$, where $`c_N`$ is the partition function of the measure $`\mu _N`$ in the case $`a_k=1,k=1,\mathrm{},N`$. Therefore, we can apply our results for $`p=1`$ to the number of summands in the random composition drawn according to $`\lambda _N`$.

Acknowledgement

The research of the second author was supported by the Fund for the Promotion of Research at Technion.

The authors are thankful to Prof. Andrew Barbour for his remarks on the first draft of the paper and to referees for their stimulating and constructive suggestions.


format: text

text:
# 1 Introduction

## 1 Introduction

The Jacobi version of the Maupertuis principle establishes that the dynamics in a natural system can be viewed as geodesic motion in an associated Riemannian manifold: If $`S=(TU)𝑑t`$ is the (natural) action functional of a system defined in a Riemannian manifold, $`(M,g)`$, the critical trajectories of $`S`$ with energy $`E=T+U`$ coincide with the extremals (geodesics) of the length functional $`L^J=𝑑s`$ defined in $`(M,h)`$, where $`h`$ is the Jacobi metric $`h=2(EU)g`$.

In this work we analyze whether or not this idea works at the level of second-order variational calculus. Namely, can the Hessian operator of a natural dynamical system be viewed as the geodesic deviation operator of the associated Riemannian manifold? Early on,it was pointed out that the answer is negative  and consequently the stability criterion provided by the geodesic deviation operator (“geometric stability criterion”) does not coincide with the usual criterion based on analysis of the spectrum of the Hessian operator (“dynamical stability criterion”). We shall show that even in this situation the Morse theory associated with the dynamical problem is the same as the Morse theory built on the geodesics of the Jacobi metric. The crux of the matter is that the Jacobi weak stability criterion provided by the number of conjugate points crossed by a trajectory is the same from both points of view.

There are several topics where the geometric criterion has been used, e.g. study of the chaotic behaviour of Hamiltonian systems  or the non-integrability of dynamical systems . The Jacobi-Maupertuis principle has been extended to the case of Lorentzian manifolds in .

## 2 The Maupertuis-Jacobi Principle

### 2.1 Preliminaries and Notation

Let $`M(M,g)`$ be a Riemannian manifold. Geodesics in $`M`$ are extremals of the free-action $`S_0`$ or length functionals $`L`$ for any differentiable curve $`\gamma :[t_1,t_2]M`$:

$$S_0[\gamma ]=_{t_1}^{t_2}\frac{1}{2}\dot{\gamma }(t)^2𝑑t;L[\gamma ]=_{t_1}^{t_2}\dot{\gamma }(t)𝑑t.$$
(1)
The critical paths are the solutions of the Euler-Lagrange equations:

$$\delta S_0=0_{\dot{\gamma }}\dot{\gamma }=0;\delta L=0_{\dot{\gamma }}\dot{\gamma }=\lambda (t)\dot{\gamma },\lambda (t)=\frac{d^2t}{ds^2}\left(\frac{ds}{dt}\right)^2.$$
(2)
As a consequence of the invariance of the length functional under reparametrizations, variational calculus on $`L`$ leads to equations where the geodesics are parametrized by an arbitrary parameter $`t`$ (often called pre-geodesics).

To decide whether a critical path is a local minimum we focus on the second variation functionals:

$$\delta ^2S_0=_{s_1}^{s_2}\mathrm{\Delta }V,V𝑑s,\delta ^2L=_{s_1}^{s_2}\mathrm{\Delta }V^{},V^{}𝑑s$$

$$\mathrm{\Delta }V=_\gamma ^{}_\gamma ^{}VR(\gamma ^{},V)\gamma ^{}=\frac{\mathrm{D}^2V}{ds^2}K_\gamma ^{}(V).$$

Here, $`V`$ denotes a proper variation ($`V(s_1)=V(s_2)=0`$, see for instance  for details) and $`V^{}`$ is the component of $`V`$ orthogonal to the geodesic. $`R`$ and $`K`$ are respectively the curvature and sectional curvature tensors. Thus, $`\mathrm{\Delta }`$ is the geodesic deviation operator.

In a natural dynamical system where the configuration space is the Riemannian manifold $`(M,g)`$, the dynamics is governed by the action:

$$S[\gamma ]=_{t_1}^{t_2}\left(\frac{1}{2}\dot{\gamma }^2U(\gamma (t))\right)𝑑t,$$
(3)
and the Euler-Lagrange-Newton equations are:

$$\delta S=0\frac{\mathrm{D}\dot{\gamma }}{dt}=\mathrm{grad}U.$$
(4)
The local stability of a given extremal under proper variations is determined by the Hessian or second-variation functional:

$$\delta ^2S=_{t_1}^{t_2}𝑑t\left(\mathrm{\Delta }V,VH(U)(V,V)\right)=_{t_1}^{t_2}𝑑t\mathrm{\Delta }V_V\mathrm{grad}U,V.$$
(5)
and we shall denote the differential operator in the quadratic form (5) as: $`\overline{\mathrm{\Delta }}V=\mathrm{\Delta }V_V\mathrm{grad}U`$.

### 2.2 Geodesics in the Jacobi Metric

Let $`x^i`$ be a system of local coordinates in $`M`$, and let us set $`ds_g^2=_{ij}g_{ij}dx^idx^j_{ij}g_{ij}dx^idx^j`$. The Jacobi metric $`h`$ is defined as a conformal transformation of $`g`$: $`ds_h^2=_{ij}h_{ij}dx^idx^j`$, $`h_{ij}=2(i_1U)g_{ij}`$, where $`i_1`$ is a constant and the Riemannian character of $`h`$ restricts the admissible values for $`i_1`$ by means of the inequality $`U<i_1`$.

Given any two vector fields $`X,Y\mathrm{\Gamma }(TM)`$, we shall use the following notation: $`h(X,Y)=X,Y^J`$, $`X^J=\sqrt{X,X^J}`$, and we shall write $`s`$ for $`s_h`$ for simplicity: $`dsds_h`$.

Theorem of Jacobi. The extremal trajectories of the variational problem associated with the functional (3) with energy $`E=i_1`$, are the geodesics of the manifold $`(M,h)`$, where $`h`$ is the Jacobi metric: $`h=2(i_1U)g`$.

The proof of this theorem can be found in several References (see  for instance). From an analytic point of view, the theorem simply establishes that the Newton equations (4) arising from the action $`S`$ are tantamount to the geodesic equations in $`(M,h)`$:$`_\gamma ^{}^J\gamma ^{}=0`$, $`\gamma ^{}=\frac{d\gamma }{ds}`$. The equivalence is shown by performing both a conformal transformation, $`h=2(i_1U)g`$, and a reparametrization from the dynamical time $`t`$ to the arc-length parameter $`s`$ in $`(M,h)`$,

$$\frac{ds}{dt}=2\sqrt{(i_1U(\gamma (s)))T}=(i_1U(\gamma (s))),$$
(6)
over the solutions (note that in the domain $`U<i_1`$ the reparametrization is well defined).

## 3 “Geometric Stability”versus “Dynamical Stability”

In the language of Variational Calculus, the Jacobi theorem states that $`\delta S=0\delta S_0^J=0`$, or, for geodesics parametrized with respect to the arc-length, $`\delta S=0\delta L^J=0`$. In  we showed the following theorems referring to the second-variation functionals:

Theorem 1. Let $`\gamma (t)`$ be an extremal of the $`S[\gamma ]`$ functional, and let $`S_0^J[\gamma ]`$ be the free-action functional for the Jacobi metric associated with $`S[\gamma ]`$. The corresponding Hessian functionals satisfy:

$$\delta ^2S_0^J[\gamma ]=\delta ^2S[\gamma ]+_{t_1}^{t_2}𝑑t\mathrm{\hspace{0.17em}2}\dot{\gamma },\frac{\mathrm{D}V}{dt}F,V,$$
(7)
where $`F=\mathrm{grad}\mathrm{ln}(2(i_1U))`$.

Theorem 2. Let $`\gamma `$ be an extremal of the $`S[\gamma ]`$ functional and let $`L^J[\gamma ]`$ be the length functional of the Jacobi metric associated with $`S[\gamma ]`$. The corresponding Hessian functionals satisfy:

$$\delta ^2L^J[\gamma ]=\delta ^2S[\gamma ]_{t_1}^{t_2}\frac{dt}{2(i_1U)}\left[_{\dot{\gamma }}\dot{\gamma },V\dot{\gamma },_{\dot{\gamma }}V\right]^2$$
(8)
Equation (8) shows that geodesics of minimal Jacobi length $`L^J[\gamma ]`$ are equivalent to stable solutions of the dynamical system, although the converse is not necessarily true. Allowing only orthogonal variations, $`V=V^{}`$, (8) can be re-written as:

$$\delta ^2S|_{V=V^{}}=\delta ^2L^J+_{s_1}^{s_2}𝑑s\left(F^J,V^{}^J\right)^2$$

The proofs of these two theorems are based on the behaviour of the covariant derivatives and the curvature tensor under reparametrizations and conformal transformations of the metric tensor.

## 4 Jacobi Fields and Morse Theory

A basic tool for the analysis of the stability properties of a given extremal of a functional is the study of the Jacobi fields: the elements of the Kernel of the second-variation or Hessian operator. Direct computation of the Jacobi fields in a given problem is in general a very complicated task: one needs to solve the Euler-Lagrange equations, calculate the second-variation operator and, finally, determine the Kernel. Fortunately, there exists a short-cut if one knows a parametric family of extremals (see, for instance,  for details):

Proposition. If $`\gamma =\gamma (t,a)`$ is a family of extremals of $`S[\gamma ]`$ characterized by the value of the real parameter $`a`$, then the vector field $`\frac{\gamma }{a}`$ is a Jacobi field.

Knowledge of the Jacobi fields allows us to calculate the conjugate points along an extremal $`\gamma `$, and it is thus possible to apply the Jacobi criterion for weak minimizers of a functional which counts the number of conjugate points crossed by the extremal as a measure of the instability.

Morse Theory establishes a link between the topology of a manifold and the critical-point structure of the functions defined in the manifold . We shall deal with Morse Theory in spaces of closed paths with a fixed point in Riemannian manifolds. The aim is the application to the separatrix trajectories in the Garnier System. A formulation of this theory à la Bott  is as follows:

Let $`\mathrm{\Omega }M=\{\gamma :S^1M/\gamma (0)=\gamma (1)=m_0\}`$ be the loop space in $`M`$ with a fixed base point. Then, some of the topological properties of $`\mathrm{\Omega }M`$ are codified in the Poincaré series $`P_t\left(\mathrm{\Omega }M^n\right)=_{k=0}^{\mathrm{}}b_kt^k`$, where $`b_k=\mathrm{dim}H_k(\mathrm{\Omega }M,)`$ are the corresponding Betti numbers.

For any functional $`S`$ defined on $`\mathrm{\Omega }M`$, we define the Morse series as $`_t\left(S\right)=_{N_c}P_t\left(N_c\right)t^{\mu (N_c)}`$, where the sum is over the critical manifolds $`N_c`$, formed either by a single isolated critical path or by a continuous set of critical paths. $`\mu (\gamma _c)`$, the Morse index, is the dimension of the subspace where the Hessian operator of the functional $`S`$, $`\mathrm{\Delta }`$, along $`\gamma _c`$ is negative definite.

The Morse inequalities, $`_t(S)P_t\left(\mathrm{\Omega }M\right)`$, tell us that the topology of $`\mathrm{\Omega }M`$ forces the existence of most of the extremals of the functional $`S`$.

The Morse index can be computed using the Morse index Theorem: The Morse index of a critical path $`\gamma _c`$ is equal to the number of conjugate points to the base point crossed by $`\gamma _c`$ counted with multiplicity. Note that this theorem is equivalent to the Jacobi criterion for weak minimizers.

From this point of view, knowledge of a one-parametric family of solutions $`\gamma (t;a)`$ of the dynamical problem (3) informs us about the geodesics $`\gamma (s;a)`$ in $`(M,h)`$. There is a one-to-one correspondence between the trajectories of the dynamical system and the geodesics of the Jacobi metric, which only differ by proper reparametrizations. Therefore, the two kinds of extremals share all the conjugate points, and the structure of the critical points of $`S`$ and $`L^J`$ is the same. The Morse series coincide, $`(S)=(L^J)`$, and the Morse theories built from the trajectories in $`(M,g)`$ and the geodesics in $`(M,h)`$ are identical.

## 5 Separatrix trajectories in the Garnier System

As an application, we shall analyze the separatrix trajectories in the Garnier System . The configuration space is $`(^N,\delta _{ij})`$ and the action functional reads:

$$S=𝑑t\left\{\frac{1}{2}\dot{\stackrel{}{q}}^2+\frac{1}{2}\left(\stackrel{}{q}\stackrel{}{q}1\right)^2+\underset{i=1}{\overset{N}{}}\frac{1}{2}\sigma _i^2q_i^2\right\}.$$
(9)
Here $`\sigma _i`$ are constants, such that $`\sigma _i\sigma _j`$, $`ij`$, and the Euler-Lagrange equations are:

$$\frac{d^2q_i}{dt^2}=2q_i\left(\stackrel{}{q}\stackrel{}{q}1\right)+\sigma _i^2q_i,i=1,2,\mathrm{},N.$$
(10)
The system is completely integrable and, in fact, Garnier found all the periodic solutions in terms of hyperelliptic functions. In  we computed the separatrix solutions lying on the boundary between bounded and unbounded motion . We shall focus in the $`N=2`$ case for simplicity (in  the very rich $`N=3`$ case is analyzed thoroughly).

For any path $`\stackrel{}{q}(t)=(q_1(t),q_2(t))`$ in $`(^2,\delta _{ij})`$, the Lagrangian $`=(q_1,q_2,`$ $`\dot{q}_1,\dot{q_2})`$ is:

$$=\frac{1}{2}\left(\dot{q}_1^2+\dot{q}_2^2\right)+\frac{1}{2}\left(q_1^2+q_2^21\right)^2+\frac{\sigma ^2}{2}q_2^2.$$
(11)
With no loss of generality, we have set $`\sigma _1=0`$, $`\sigma _2=\sigma `$. We shall only deal with the case $`\sigma <1`$, because the separatrix solutions in this regime have a richer structure. In fact $`\sigma >1`$ generates a situation where the separatrices are reduced to a single singular solution. The system is not only completely integrable, but Hamilton-Jacobi separable in Jacobi elliptical coordinates: $`q_1^2=\sigma ^2(1\mu _1)(1\mu _2)`$, $`q_2^2=\sigma ^2(\overline{\sigma }^2\mu _1)(\overline{\sigma }^2\mu _2)`$, where $`\overline{\sigma }^2=1\sigma ^2`$. With this change, $`^2`$ is mapped into the infinite parallelogram determined by the inequalities: $`\mathrm{}<\mu _1<\overline{\sigma }^2<\mu _2<1`$.

The choice of the separatrix trajectories is achieved by setting the values of the two invariants in involution to zero: $`i_1=i_2=0`$. $`i_1`$ denotes the energy first-integral and $`i_2`$ is the second first-integral of the system, functionally independent of $`i_1`$ (see for instance  for details). This choice confines the dynamics to the finite parallelogram $`\overline{𝐏}(0)0\mu _1<\overline{\sigma }^2<\mu _2<1`$ (in Cartesian coordinates it is the area bounded by the ellipse $`q_1^2+\frac{q_2^2}{\overline{\sigma }^2}=1`$).

In elliptical coordinates the Euclidean metric $`\delta _{ij}`$ becomes : $`g_{12}=0=g_{21}`$, $`g_{11}=\frac{(\mu _1\mu _2)}{4(\mu _11)(\mu _1\overline{\sigma }^2)}`$, $`g_{22}=\frac{(\mu _2\mu _1)}{4(\mu _21)(\mu _2\overline{\sigma }^2)}`$. The Christoffel symbols are obtained from

$`\mathrm{\Gamma }_{11}^1={\displaystyle \frac{1}{2(\mu _1\mu _2)}}{\displaystyle \frac{2\mu _1(1+\overline{\sigma }^2)}{2(\mu _11)(\mu _1\overline{\sigma }^2)}};\mathrm{\Gamma }_{12}^1={\displaystyle \frac{1}{2(\mu _1\mu _2)}};\mathrm{\Gamma }_{22}^1={\displaystyle \frac{(\mu _11)(\mu _1\overline{\sigma }^2)}{2(\mu _21)(\mu _2\overline{\sigma }^2)(\mu _1\mu _2)}}`$
Replacing 1 by 2, $`\mathrm{\Gamma }_{22}^2`$ is equivalent to $`\mathrm{\Gamma }_{11}^1`$, and the same happens between $`\mathrm{\Gamma }_{12}^2`$, $`\mathrm{\Gamma }_{11}^2`$, and $`\mathrm{\Gamma }_{12}^1`$, $`\mathrm{\Gamma }_{22}^1`$. The rest of the symbols are trivially related to these by symmetries or are null.

The Jacobi metric associated with the system, with $`i_1=0`$, is $`h=\frac{1}{\mu _2\mu _1}((\mu _1^3\overline{\sigma }^2\mu _1^2)+(\mu _2^3\overline{\sigma }^2\mu _2^2))g`$, and the Christoffel symbols can be obtained by direct calculation ().

### 5.1 Singular Solutions

There are two singular solutions to the system, and these live on the border of the parallelogram $`\overline{𝐏}(0)`$.

Taking $`q_2=0`$, the equations of the motion can be easily integrated, even in Cartesian coordinates, to find the solutions: $`\stackrel{}{q}(t)=(\overline{q}_1(t),0)=(\pm \mathrm{tanh}(tt_0),0)`$. In elliptical coordinates, these solutions lie on two edges of the border of $`\overline{𝐏}(0)`$: I, $`\stackrel{}{\mu }_\mathrm{I}(t)=(1\mathrm{tanh}^2(\pm t),\overline{\sigma }^2)`$, for $`t(\mathrm{},arctanh(\sigma )][arctanh(\pm \sigma ),\mathrm{})`$, and II, $`\stackrel{}{\mu }_{\mathrm{II}}(t)=(\overline{\sigma }^2,1\mathrm{tanh}^2(t))`$ for $`t[arctanh(\sigma ),`$ $`arctanh(\pm \sigma )]`$ (taking $`t_0=0`$).

The geodesic equations of the Jacobi metric reduced to the $`q_2=0`$ orbit become a single non-trivial differential equation, which is solved by a cubic algebraic equation :

$$q_1^{\prime \prime }+\frac{2q_1}{q_1^21}q_1^2=0\overline{q}_1\frac{\overline{q}_1^3}{3}=\pm s.$$
(12)
$`(t_1,t_2)=(\mathrm{},\mathrm{})`$ has been reparametrized to $`[s_1,s_2]=[\frac{2}{3},\frac{2}{3}]`$. The explicit solution of the cubic (12) is

$$q_1(s)=\overline{q}_1(s)=\mathrm{cos}\frac{\theta }{3}+\sqrt{3}\mathrm{sin}\frac{\theta }{3},$$
(13)
where $`\theta =\mathrm{arctan}\frac{\sqrt{49s^2}}{3s}`$ .

There exists a second type of singular solutions living on the edge $`\mu _1=0`$, which also belongs to the border of $`\overline{𝐏}(0)`$. This edge becomes the ellipse $`q_1^2+\frac{q_2^2}{\overline{\sigma }^2}=1`$ in Cartesian coordinates. The solutions are either $`\stackrel{}{q}(t)=(\mathrm{tanh}(\pm \sigma t),\overline{\sigma }sech(\pm \sigma t))`$ in $`^2`$, or $`\stackrel{}{\mu }(t)=(0,1\sigma ^2\mathrm{tanh}^2\sigma t)`$ in the elliptic plane.

The geodesic equations of the Jacobi metric reduced to the $`q_1^2+\frac{q_2^2}{\overline{\sigma }^2}=1`$ orbit also become a single non-trivial differential equation, which is solved by a cubic algebraic equation :

$$q_1^{\prime \prime }\frac{2\sigma ^2q_1}{1\sigma ^2q_1^2}q_1^2=0\sigma \left(\overline{q}_1\frac{\sigma ^2\overline{q}_1^3}{3}\right)=s.$$
(14)
The explicit solution of the cubic together with the elliptic orbit afford the geodesic:

$`\overline{q}_1(s)={\displaystyle \frac{1}{\sigma }}\mathrm{cos}{\displaystyle \frac{\theta }{3}}+{\displaystyle \frac{\sqrt{3}}{\sigma }}\mathrm{sin}{\displaystyle \frac{\theta }{3}},\overline{q}_2(s)={\displaystyle \frac{\overline{\sigma }}{\sigma }}\sqrt{2+\sigma ^2+\mathrm{cos}{\displaystyle \frac{2\theta }{3}}+\sqrt{3}\mathrm{sin}{\displaystyle \frac{2\theta }{3}}}`$
from $`s_1=\sigma (1\frac{\sigma ^2}{3})`$ to $`s_2=\sigma (1\frac{\sigma ^2}{3})`$.

### 5.2 General Solution

The Hamiltonian associated with the Lagrangian (11), in elliptical coordinates is of the Stäckel type:

$`H=`$ $`{\displaystyle \frac{1}{2(\mu _1\mu _2)}}(4(\mu _11)(\mu _1\overline{\sigma }^2)\pi _1^24(1\mu _2)(\mu _2\overline{\sigma }^2)\pi _2^2`$ (15)
$`(\mu _1^3\overline{\sigma }^2\mu _1^2)+(\mu _2^3\overline{\sigma }^2\mu _2^2))`$
and the integration of the Hamilton-Jacobi equation, for separation constants equal to zero, , leads to the equation of the orbits,

$`\left(\right|{\displaystyle \frac{\sqrt{1\mu _1}\sigma }{\sqrt{1\mu _1}+\sigma }}||{\displaystyle \frac{\sqrt{1\mu _1}+1}{\sqrt{1\mu _1}1}}|^\sigma )^{\mathrm{sign}(\pi _1)}`$
$`\left(\left|{\displaystyle \frac{\sqrt{1\mu _2}\sigma }{\sqrt{1\mu _2}+\sigma }}\right|\left|{\displaystyle \frac{\sqrt{1\mu _2}+1}{\sqrt{1\mu _2}1}}\right|^\sigma \right)^{\mathrm{sign}(\pi _2)}=e^{2\sigma \overline{\sigma }^2a},`$ (16)
and the equation of the temporal dependence,

$$\left|\frac{\sqrt{1\mu _1}\sigma }{\sqrt{1\mu _1}+\sigma }\right|^{\overline{\sigma }^2\mathrm{sign}(\pi _1)}\left|\frac{\sqrt{1\mu _2}\sigma }{\sqrt{1\mu _2}+\sigma }\right|^{\overline{\sigma }^2\mathrm{sign}(\pi _2)}=e^{2(t+t_0)\sigma }.$$
(17)
Here $`a`$ is a constant that parametrizes the different orbits and $`t_0`$ is another constant coming from the time-translation invariance.

At the limit $`a\pm \mathrm{}`$, equation (16) reduces to the singular solution $`q_2=0`$ together with the singular solution living on the“ellipse”. Thus, identifying $`a=\mathrm{}`$ with $`a=\mathrm{}`$, we see that all the orbits are parametrized by a periodic parameter $`a`$ : the space of all orbits is the $`S^1`$ circle.

A straightforward application of the Hamilton-Jacobi procedure to the free-Hamiltonian, $`H^J=\frac{1}{2}\left(h^{11}\stackrel{~}{\pi }_1^2+h^{22}\stackrel{~}{\pi }_2^2\right)`$ , $`\stackrel{~}{\pi }_j=\frac{d\mu _j}{ds}`$, provides all the geodesics of the Jacobi metric with $`i_1=0`$.

The equation giving the orbits is exactly (16), and the arc-length dependence is determined by the expression:

$$s+s_0=(1)^{\mathrm{sign}(\pi _1)}\frac{\mu _1+2}{3}\sqrt{1\mu _1}(1)^{\mathrm{sign}(\pi _2)}\frac{\mu _2+2}{3}\sqrt{1\mu _2}$$
(18)
### 5.3 Jacobi fields

In  the Jacobi fields are calculated on the trajectories. Here, we closely follow the same idea to compute the Jacobi fields on the geodesics, which are the solutions of the equations (16)-(18). Thus, we have a two-parametric family $`\stackrel{}{\mu }(s;a,s_0)`$ of geodesics, and, accordingly, two Jacobi fields: $`\frac{\stackrel{}{\mu }}{s_0}`$ and $`\frac{\stackrel{}{\mu }}{a}`$; see Section §2.2 . The first Jacobi field is tangent to the geodesic because it obeys the invariance under $`s`$ translations. We focus on the computation of the orthogonal Jacobi field: $`J=\frac{\stackrel{}{\mu }}{a}`$. An explicit formula for the geodesic $`\stackrel{}{\mu }(s;a,s_0)=(\overline{\mu }_1,\overline{\mu }_2)`$ in terms of known analytical functions is not available, but we implicitly derive the system (16)-(18) with respect to $`a`$ to obtain a linear system in the components of $`J`$:

$`{\displaystyle \frac{(1)^{\alpha _1}}{\overline{\mu }_1\sqrt{1\overline{\mu }_1}(\overline{\sigma }^2\overline{\mu }_1)}}{\displaystyle \frac{\overline{\mu }_1}{a}}{\displaystyle \frac{(1)^{\alpha _2}}{\overline{\mu }_2\sqrt{1\overline{\mu }_2}(\overline{\sigma }^2\overline{\mu }_2)}}{\displaystyle \frac{\overline{\mu }_2}{a}}=2`$ (19)
$`{\displaystyle \frac{(1)^{\alpha _1}\overline{\mu _1}}{\sqrt{1\overline{\mu }_1}}}{\displaystyle \frac{\overline{\mu }_1}{a}}+{\displaystyle \frac{(1)^{\alpha _2}\overline{\mu }_2}{\sqrt{1\overline{\mu }_2}}}{\displaystyle \frac{\overline{\mu }_2}{a}}=0`$ (20)
The Jacobi field that solves (19)-(20) is:

$$J=j(\overline{\mu }_1,\overline{\mu }_2)\left((1)^{\alpha _1}\overline{\mu }_2\sqrt{1\overline{\mu }_1}\frac{}{\mu _1}+(1)^{\alpha _2}\overline{\mu }_1\sqrt{1\overline{\mu }_2}\frac{}{\mu _2}\right),$$
(21)
where

$$j(\overline{\mu }_1,\overline{\mu }_2)=\frac{2\overline{\mu }_1\overline{\mu }_2(\overline{\mu }_1\overline{\sigma }^2)(\overline{\mu }_2\overline{\sigma }^2)}{(\overline{\mu }_1\overline{\mu }_2)(\overline{\mu }_1^2+\overline{\mu }_2^2+\overline{\mu }_1\overline{\mu }_2(\overline{\mu }_1+\overline{\mu }_2)\overline{\sigma }^2)},\alpha _1,\alpha _2=0,1$$

$`J`$ is zero both at the starting point D and at the focus F of the ellipse, and, henceforth, the focus is a conjugate point of D of multiplicity 1.

### 5.4 The Morse series

Taking $`M=\overline{𝐏}_2(0)`$; $`\mathrm{\Omega }M=𝒞`$, we have the following critical point structure for the geodesics that start and end at the point D. Firstly, the point D itself is a possible geodesic with a Morse index 0. Secondly, all members of the family determined by the general solution, $`\stackrel{}{\mu }(s;a,s_0)`$ cross the focus. Bearing in mind that $`P_t(S^1)=(1+t)`$, the family contributes to the Morse series with $`(1+t)t`$. Iterating two of the solutions, $`\stackrel{}{\mu }\mathrm{}\stackrel{}{\mu }`$, D is a conjugate point to itself and the contribution to the Morse series is $`(1+t)t^3`$. More iterations tell us that $`_t(L^J|_𝒞)=1+(1+t)t+(1+t)t^3+\mathrm{}=\frac{1}{1t}`$. The Morse series is equal to the Poincaré series of the space of closed geodesics in $`S^2`$: $`_t\left(S|_𝒞\right)=P_t\left(\mathrm{\Omega }S^2\right)`$.

In  we calculated the Morse series for the $`N=3`$ case. The result is: $`_t\left(L^J|_𝒞\right)=P_t\left(𝒞\right)=\frac{1}{1t^2}=P_t\left(\mathrm{\Omega }S^3\right)`$.


text:
# Incompressibility of strange matter ∗.

## 1 Introduction

An exciting issue of modern astrophysics is the possible existence of a family of compact stars made entirely of deconfined u,d,s quark matter or “strange matter” (SM) and thereby denominated strange stars (SS). They differ from neutron stars, where quarks are confined within neutrons, protons and eventually within other hadrons (hadronic matter stars). The possible existence of SS is a direct consequence of the so called strange matter hypothesis, according to which the energy per baryon of SM would be less than the lowest energy per baryon found in nuclei, which is about 930 MeV for $`Fe^{56}`$. Also, the ordinary state of matter, in which quarks are confined within hadrons, is a metastable state. Of course, the hypothesis does not conflict with the existence of atomic nuclei as conglomerates of nucleons, or with the stability of ordinary matter.

The best observational evidence for the existence of quark stars come from some compact objects, the X-Ray burst sources SAX J1808.4$``$3658 (the SAX in short) and 4U 1728$``$34, the X-ray pulsar Her X-1 and the superburster 4U 1820$``$30. The first is the most stable pulsating X-ray source known to man as of now. This star is claimed to have a mass $`M1.4M_{}`$ and a radius of about 7 kms . Coupled to this claim are the various other evidences for the existence of SS, such as the possible explanation of the two kHz quasi-periodic oscillations in 4U 1728 - 34  and the quark-nova explanation for $`\gamma `$ ray bursts .

The expected behaviour of SS is directly opposite to that of a neutron star as Fig.(1) shows. The mass of 4U 1728$``$34 is claimed to be less than 1.1 $`M_{}`$ in Li et al. , which places it much lower in the M-R plot and thus it could be still gaining mass and is not expected to be as stable as the SAX. So for example, there is a clear answer  to the question posed by Franco : why are the pulsations of SAX not attenuated, as they are in 4U 1728$``$34 ?

From a basic point of view the equation of state for SM should be calculated solving QCD at finite density. As we know, such a fundamental approach is presently not feasible even if one takes recourse to the large colour philosophy of ’t Hooft . A way out was found by Witten  when he suggested that one can borrow a phenomenological potential from the meson sector and use it for baryonic matter. Therefore, one has to rely on phenomenological models. In this work, we use different equations of state (EOS) of SM proposed by Dey et al  using the phenomenological Richardson potential. Other variants are now being proposed, for example the chromo dielectric model calculations of Malheiro et al..

Fig.(2) shows the energy per baryon for the EOS of . One of them (eos1, SS1 of ) has a minimum at E/A = 888.8 $`MeV`$ compared to 930.4 of $`Fe^{56}`$, i.e., as much as 40 $`MeV`$ below. The other two have this minimum at 911 $`MeV`$ and 926 $`MeV`$, respectively, both less than the normal density of nuclear matter. The pressure at this point is zero and this marks the surface of the star in the implementation of the well known TOV equation. These curves clearly show that the system can fluctuate about this minimum, so that the zero pressure point can vary.

## 2 Incompressibility : its implication for Witten’s Cosmic Separation of Phase scenario.

In nuclear physics incompressibility is defined as

$$K=9\frac{}{n}\left(n^2\frac{\epsilon }{n}\right),$$
(1)
where $`\epsilon =E/A`$ is the energy per particle of the nuclear matter and $`n`$ is the number density. The relation of K with bulk modulus B is

$$K=\frac{9B}{n}.$$
(2)
$`K`$ has been calculated in many models. In particular, Bhaduri et al  used the non - relativistic constituent quark model, as well as the bag model, to calculate $`K_A`$ as a function of $`n`$ for the nucleon and the delta. They found that the nucleon has an incompressibility $`K_N`$ of about 1200 MeV, about six times that of nuclear matter. They also suggested that at high density $`K_A`$ matches onto quark gas incompressibility.

The velocity of sound in units of light velocity c is given by

$$v=\sqrt{K/9\epsilon }.$$
(3)
The simple models of quark matter considered in  use a Hamiltonian with an interquark potential with two parts, a scalar component (the density dependent mass term) and a vector potential originating from gluon exchanges. In the absence of an exact evaluation from QCD, this vector part is borrowed from meson phenomenology . In common with the phenomenological bag model, it has built in asymptotic freedom and quark confinement (linear). In order to restore the approximate chiral symmetry of QCD at high densities, an ansatz is used for the constituent masses, viz.,

$$M(n)=M_Qsech\left(\nu \frac{n}{n_0}\right),$$
(4)
where $`n_0`$ is normal nuclear matter density and $`\nu `$ is a parameter. There may be several EOS’s for different choices of parameters employed to obtain the EOS. Some of them are given in the table (1). There $`\alpha _s`$ is perturbative quark gluon coupling and $`\mathrm{\Lambda }`$ is the scale parameter appeared in the vector potential. Changing other parameters one can obtain more EOS’s. The table (1) also shows the masses($`M_G`$),radii (R) and the number density($`n_s`$) at star surface of maximum mass strange stars obtained from the corresponding EOS’s. The surface of the stars occur at the minimum $`\epsilon `$ where pressure is zero.

The general behaviour of the curves is relatively insensitive to the parameter $`\nu `$ in $`M(n)`$ as well as the gluon mass, as evident from the figure(3).

In figure (4), $`K`$ with three values of $`M_Q`$ implying different running masses, $`M(n)`$, is plotted as a function of the density expressed by its ratio to $`n_0`$. Given for comparison, is the incompressibility $`K_q`$ of a perturbative massless three flavour quark gas consisting of zero mass current quarks  using the energy expression given in  to order $`\alpha _s^2`$.

It can be seen that as $`M_Q`$ decreases, the nature of the relation approaches the perturbative case of . At high density our incompressibility and that due to Baym  matches, showing the onset of chiral symmetry restoration. In EOS1 for uds matter, the minimum of $`\epsilon `$ occurs at about $`4.586n_0`$. nucleation may occur at a density less than this value of $`n`$. This corresponds to a radius of about $`0.67fm`$ for a baryon. For EOS1 we find $`K`$ to be $`1.293GeV`$ per quark at the star surface.

It is encouraging to see that this roughly matches with the compressibility $`K_N`$ so that no ‘phase expands explosively’. In the Cosmic Separation of Phase scenario, Witten  had indicated at the outset that he had assumed the process of phase transition to occur smoothly without important departure from equilibrium. If the two phases were compressed with significantly different rates, there would be inhomogenieties set up.

But near the star surface at $`n4`$ to $`5n_0`$ the matter is more incompressible showing a stiffer surface. This is in keeping with the stability of strange stars observed analytically with the Vaidya-Tikekar metric by Sharma, Mukherjee, Dey and Dey .

The velocity of sound, $`v_s`$ peaks somewhere around the middle of the star and then falls off. We show it for three different EOS in fig.(5) with parameters given in Table (1).

Next we turn to the model of Zimanyi and Moskowski , where the coupling constants are density dependent. It was shown  that the quark condensate derived from this model via the Hellmann-Feynman theorem is physically acceptable in this model. Details may be found in  but the essence is that in the more conventional Walecka model the condensate increases with density whereas in QCD a decrease is expected.

In a recent paper  Krein and Vizcarra (KV in short) have put forward an EOS for nuclear matter which exhibits a transition from hadronic to quark matter. KV start from a microscopic quark-meson coupling Hamiltonian with a density dependent quark-quark interaction and construct an effective quark-hadron Hamiltonian which contains interactions that lead to quark deconfinement at sufficiently high densities. At low densities, their model is equivalent to a nuclear matter with confined quarks, i.e., a system of non-overlapping baryons interacting through effective scalar and vector meson degrees of freedom, while at very high densities it is ud quark matter. The $`K_{NM}`$ at the saturation density is fitted to be $`248MeV`$. This EOS also gives a smooth phase transition of quark into nuclear matter and thus, conforms to Witten’s assumption. Interestingly enough, the transition takes place at about $`5n_0`$.

The KV model does not incorporate strange quarks so that comparison with our EOS is not directly meaningful. However it is quite possible that the signal results of the KV calculations mean that quark degrees of freedom lower the energy already at the ud level and once the possibility of strange quarks is considered the binding exceeds that of $`Fe^{56}`$. An extension of KV with strange quarks is in progress111G. Krein, by e-mail.

Results obtained from the KV calculation are presented simultaneously. The incompressibility shows softening and the velocity of sound decreases when quark degrees of freedom open up, as expected. At $`5n_0`$ $`K_{qN}`$ is about $`2GeV`$ in qualitative agreement with our value.

Note that in our EOS we also have strange quarks reducing the value of $`K`$. Fig (9) shows the $`v_s`$ as a function of $`n_B/n_0`$ for both kinds of EOS. The EOS with quarks shows a lowering of $`v_s`$.

## 3 Discussions.

In our model the density is large even at the surface ($`4.586n_0`$ for $`EOS1`$) so that the quarks in the fermi sea have high momentum on the average. Therefore, they often come close together. It was shown by Bailin and Love that in certain colour channels , the short range one gluon exchange spin-spin interaction can produce extra attraction, forming diquarks. In s-state the pair can be formed only in flavour antisymmetric state. The flavour antisymmetric diquark may form in colour symmetric(6) and spin symmetric (triplet) state, or in colour antisymmetric ($`\overline{3}`$) and spin antisymmetric (singlet) state. The potential strength of the latter is 6 times the former. The N-$`\mathrm{\Delta }`$ mass splitting is obtainable from the one gluon exchange potential, the delta function part of which is smeared to a Gaussian,

$$H_{i,j}=80.51\sigma ^3(\lambda _i.\lambda _j)(S_i.S_j)e^{\sigma ^2r_{ij}^2}.$$
(5)
In the above, the strength of this spin dependent part is taken from Dey and Dey(1984) . The pairing energy of ud pairs in the spin singlet and colour $`\overline{3}`$ state is $`3.84MeV`$  222With this pairing energy superbursts observed from some astrophysical x-ray sources has been explained very well.. Our model does not predict diquarks that are permanent - since the quarks must have high momentum transfer if they are to interact strongly with a force that is short-range. The formation and breaking of pairs give rise to endothermic and exothermic processes leading to fast cooling and superbursts respectively.

We also note that the lowering of energy is not very large, as compared to the approximately 42 $`MeV`$ energy difference per baryon (930.6 $`Fe^{56}`$ to $`888MeV`$), seen in EOS1. So we do not expect any drastic change in the incompressibility or the velocity of sound due to diquark pairing, nor do we expect a phase transition to a colour-flavour locked state.

We are grateful to the referee for inviting a comment on quark pairing in the strange matter scenario.


text:
# Novel aspects and strong correlation in the electronic structure of Sr2FeMoO6

## Abstract

We investigate the electronic structure of Sr2FeMoO6 combining photoemission spectroscopy with a wide range of photon energies and electronic structure calculations based on first-principle as well as model Hamiltonian approaches to reveal several interesting aspects. We find evidence for unusually strong Coulomb correlation effects both in the Fe 3$`d`$ and O 2$`p`$ states, with an enhanced manifestation in the majority spin channel. Additionally, O 2$`p`$ states exhibit a spin-splitting of nonmagnetic origin, which nevertheless is likely to have subtle influence on the stability of novel ferromagnetism of this compound.

PACS number(s): 79.60.-i, 71.20.-b, 71.27.+a, 71.20.Ps

Introduction

Sr2FeMoO6 has been recently shown to exhibit a pronounced negative magnetoresistance at higher temperatures and lower magnetic fields  compared to the manganites, making it intrinsically more suitable for technological applications. An essential aspect of understanding the properties of this material is to understand the underlying electronic structure. The basic ingredient for crystal structure of this material is the cubic perovskite structure, with Fe and Mo alternating in the (100), (010) and (001) directions. As in the case of the other $`AB`$O3 perovskites, where $`A`$ is an alkali-metal atom and $`B`$ is a transition metal atom, the valence and conduction bands of Sr2FeMoO6 involve the $`d`$ states of the two transition metal atoms - Fe and Mo, besides the oxygen $`p`$ states. It has been recently shown  that the Mo 4$`d`$-O 2$`p`$ hybridized states play a crucial role in establishing the unusually high magnetic ordering temperature in Sr2FeMoO6 via a new mechanism that relies on the location of these Mo 4$`d`$-O 2$`p`$ hybridized states energetically between the exchange-split Fe states. There has been only one direct investigation  of the electronic structure of this compound so far, with another investigation  reporting on a related compound, Ba2FeMoO6. These experimental results suggest substantial disagreement with ab initio band structure results and are interpreted in terms of LDA+U calculations by incorporating Coulomb interactions within the Fe 3$`d`$ manifold. We have performed photoemission experiments on Sr2FeMoO6 and Sr2FeMoxW1-xO6 with a photon energy range that encompasses the combined range of the previous two experiments, farther extended by including photoemission results with x-rays. These detailed experimental results are analyzed in conjunction with calculations based on multiband Hubbard model within mean-field decoupling of the interaction terms as well as ab initio approaches.

We identify the origin of various features in the experimental valence band spectrum with the help of the dependence of the matrix elements on the photon energy. Thus, we find that while the specific ordering of energy levels near the Fermi energy, crucial for the magnetism, is in agreement with the ab initio results, there are some significant differences between the experimental and the ab initio results in agreement with previous reports, indicating the presence of strong correlation effects. Moreover, we show that the conventional LDA+U and parametrized model Hamiltonian approaches are also in disagreement with the experimental results in certain details, particularly in the spectral range dominated by oxygen $`p`$ states, though the inclusion of $`U_{dd}`$ improves the description of Fe $`d`$ related spectral features as also observed previously. It turns out that the inclusion of the interaction terms within the oxygen $`p`$ states improves the agreement between the calculated results and the experiment. We further establish a few unusual aspects of the electronic structure in this compound specifically related to the O $`p`$ band, such as a spin-splitting of nonmagnetic origin. These findings are possibly of importance in understanding the stability of the ferromagnetic state and provide a clue why the magnetic structure in this compound has been recently found to depend critically on factors that influence the O $`p`$ states.

Experimental

Polycrystalline Sr2FeMoO6 was prepared by standard solid state route as described in ref. 1. Stoichiometric amounts of SrCO3, MoO3 and Fe2O3 were mixed thoroughly and heated at 900C in air for 3 hrs. The fully oxidized black powder is then grounded and pelletized into 10 mm. diameter pellets of 2 mm. thickness. These pellets were then annealed at 1200C in a flowing gas mixture of $``$10% H2/Ar for 5 hrs. X-ray diffraction experiments on this powder confirmed the presence of highly ordered (i.e. near perfect alternate occupancy of Fe and Mo ions along the three cubic axes) and pure double perovskite Sr2FeMoO6 phase. It has tetragonal I4/mmm space group with $`a`$ = $`b`$ = 10.5229 a.u. and $`c`$ = 14.9236 a.u.; Sr, Fe, Mo, O1 and O2 sites occupy 8g, 4e, 4e, 4e and 8h positions, respectively.

A different method was applied for synthesizing W-doped Sr2FeMoO6 samples, because it was otherwise not possible to entirely get rid of SrMoO4/SrWO4 type oxidized impurities, whenever the W-doped samples were synthesized by the above mentioned procedure. Therefore, all the W-doped samples were first prepared by melt-quenching method ; the technique which always produced slightly less ordered but highly pure samples. These pellets were crushed, grounded thoroughly and finally annealed at the same condition in the form of pellets in order to achieve maximum ordering.  The purity and chemical homogeneity of the final products were checked by XRD and Energy Dispersive Analysis of X-rays (EDAX) techniques.

Majority of the photoemission experiments were carried out at the VUV photoemission beamline, 3.2R at Synchrotron Radiation Centre, Elettra, Trieste, with a sample temperature of 77 K. The sample surface was cleaned by in situ scraping with a diamond file. XPS valence band spectrum was recorded in a standard VSW electron spectrometer  at $``$120 K with AlK$`\alpha `$ radiation.

Calculational Methods

The ab-initio band structure of Sr2FeMoO6 was computed using both the LMTO-ASA method using the Perdew Wang GGA exchange functional  and the plane wave pseudopotential method. No empty spheres were introduced, and the volume filling criterion was satisfied with less than 16% overlap between the spheres in the LMTO-ASA calculations. Ultrasoft pseudopotentials  as implemented in VASP  were used for the calculations performed with a plane wave cutoff of 20 Ry and the accuracy of our conclusions was checked by increasing the plane wave cutoff to 25 Ry. Both the internal as well as the external cell parameters were optimized. The GGA  optimized lattice constants were 10.59 and 14.997 a.u. and the optimized internal coordinates were Sr (8g) z=0.25, Fe (4e) z=0, Mo (4e) z=0.5, O1 (4e) z=0.2527 and O2 (8h) x=0.2474, in agreement with the experimental parameters.

The band structure obtained from both ab-initio methods were found to be in agreement. We fit the magnetic band structure along various symmetry directions with band dispersions within a nearest neighbor tight binding model  in order to obtain estimates of the hopping parameters for the use in the multiband Hubbard model. The basis of the tight-binding Hamiltonian includes $`d`$ orbitals on Fe and Mo atoms as well as $`p`$ orbitals on oxygen. As the Sr states do not contribute in the energy window that we are interested in, the basis of the tight binding model does not include states on the Sr atoms. The tight binding model considered includes hopping between oxygen atoms, as well as between Fe and oxygen and Mo and oxygen; these interaction strengths  are given in terms of the Slater Koster parametrization. Instead of using different onsite energies ($`ϵ_i`$) for up and down spin orbitals to describe the magnetic structure as in the past,  we replace the different on-site energies for the up and down spin of Fe by an on-site direct Coulomb ($`U_{dd}`$) interaction strength of 3 eV and an exchange Hund’s coupling ($`J`$) strength of 0.7 eV on the Fe site as appropriate  for Fe3+ oxide systems; this method allows us to solve for the changes in the energies of the Fe 3$`d`$ level as a result of direct and exchange Coulomb interactions self-consistently. Besides the usual one-body tight binding part, the Hamiltonian includes the intra-site Coulomb interaction term at the Fe site in the form $`(U_{dd}`$-$`J`$$`\delta _{\sigma _1,\sigma _2}`$)$`n_{\alpha _1,\sigma _1}`$ $`n_{\alpha _2,\sigma _2}`$, where $`\alpha _i`$ denotes the five 3$`d`$ orbitals on the Fe site and $`\sigma _i`$ denotes the spin. As discussed later in the text, certain features in the experimental spectra prompted us to include an additional term representating Coulomb interactions within the oxygen $`p`$ orbitals in the form $`U_{pp}`$ $`n_{\beta _1,\sigma }`$ $`n_{\beta _2,\sigma ^{^{}}}`$. Effect of such a term involving $`U_{pp}`$ has been studied in the context of the electronic structure of the of high temperature superconductors  and other oxide systems.  The four fermion terms in the Hamiltonian were decoupled by using a mean-field approximation.  The decoupled Hamiltonian involved the occupancies of various orbitals, which were then solved selfconsistently over a k-mesh of 1000 k-points.

Results and discussion

In Fig. 1 we show a selection of the valence band photoemission spectra of Sr2FeMoO6, obtained using different incident photon energies. All the spectra exhibit a finite intensity at $`E_F`$, in agreement with the metallic behavior in this compound. For comparison, we also show the calculated total density of states (DOS) for fully ordered Sr2FeMoO6, obtained within spin-polarized LMTO-ASA band structure calculation with generalized gradient approximation (GGA) over the relevant energy range just below the experimental spectra. The experimental spectra show four distinct features, appearing at about 0.4, 1.8, 5.6 and 8.3 eV binding energies, marked by the vertical dotted lines in Fig. 1. The calculated LMTO-DOS clearly does not agree with the measured spectra, with the calculation suggesting regions of low DOS at three of these energies, namely at 1.8, 5.6 and 8.3 eV. In particular, the remarkably intense feature at 8.3 eV in the experimental spectra lies outside the calculated bandwidth. This is similar to the cases of U-intermetallics  and transition metal oxides, such as NiO , exhibiting correlation-driven large spectral weights, or even distinct peaks, outside the calculated bandwidth. It is significant to note that this correlation driven satellite in Sr2FeMoO6 is very intense and is undoubtedly the highest among all Fe3+ oxides studied so far ; for example, it is known  that spectroscopic data on LaFeO3, a closely related perovskite oxide with Fe3+, can be reasonably well described within ab initio approaches.

A standard route to take into account the effect of Coulomb interactions on the band structure is to incorporate this interaction within the transition metal 3$`d`$ manifold in a mean-field sense, as in the case of LDA+U approach or the Hartree-Fock (HF) type mean-field approach described for the multiband Hubbard model in the previous section. It is indeed true that such approaches have proved to be very useful for describing transition metal compounds. We show the total density of states obtained within the HF-treatment of the multiband Hamiltonian including Coulomb interactions only at the Fe sites at the bottom of Fig. 1; our result is in agreement with LDA+U results . A comparison of the experimental data with the two calculations (LMTO and HF)readily suggests that the HF-type approach with the multiband model provides a better agreement with the experiment in comparison to the LMTO-GGA results. Specifically, the experimentally observed 1.8 eV feature and the satellite structure at about 8 eV, outside the GGA derived bandwidth, have corresponding features in the model Hamiltonian results. Interestingly, however, the sharp disagreement between the most intense feature in the 5-6 eV range in the experimental spectrum and the deep valley of low calculated DOS in the same energy range continues to exist in the HF-DOS, just as in the LMTO-DOS. The same discrepency can also be observed in the LDA+U results compared with experimental results, reported in ref.  though this disagreement has not been commented up on so far.

Before attempting to understand the origin of this discrepency, we first establish the origin of various experimental features from the dependencies of the normalized spectral intensities on the photon energies (Figs. 2(a) and (b)). We found that the two features appearing at 1.8 and 8.3 eV exhibit significant non-monotonic $`h\nu `$ dependence in the range of 42-65 eV (Fig. 2(a)), spanning the Fe 3$`p`$ absorption threshold near 55 eV. Such a resonance phenomenon clearly establishes substantial Fe 3$`d`$ contributions in these states. More interestingly, the nature of the dependencies is different for the two features. The 8.3 eV feature exhibits a clear maximum, indicating about 30% resonant enhancement, characteristic of such correlation driven satellite or incoherent features. We have found almost identical variations in the satellite features in closely related W-substituted compounds also, as shown in the figure. In contrast, the 1.8 eV feature shows approximately a 20% dip in the intensity near the on-resonant condition, as is characteristic of more band-like or coherent features.

In order to identify the states with dominant Mo contributions, we plot the normalized intensities of the features at the $`E_F`$ and 8.3 eV for $`h\nu `$ between 75 and 152 eV in Fig. 2(b). This photon energy range is most suitable for the purpose, since there is a minimum (called Cooper minimum) in the cross-section of localized Mo 4$`d`$ related intensities in solids at about $`h\nu `$ 100-110 eV.  Fig. 2(b) clearly shows a significant minimum in the relative intensity of the 8.3 eV feature (triangles) at about $`h\nu `$ 100 eV, exactly where the Cooper minimum for the localized Mo 4$`d`$ states are expected. This establishes that the Mo 4$`d`$ states also contribute significantly into the 8 eV spectral region. In contrast to the behavior of the localized states, delocalized band-like Mo $`d`$ states, however, do not exhibit any Cooper minimum.  Therefore, the relative intensity of Mo 4$`d`$ delocalized states compared to the O 2$`p`$ states is expected to show a rapid, monotonic increase with increasing photon energy, $`h\nu `$, due to the changes in the relative Mo 4$`d`$ : O 2$`p`$ cross-sections.  Such an increasing trend is indeed exhibited by the relative intensity of the feature at $`E_F`$ (see square data in Fig. 2(b)), establishing a dominant contribution from delocalized Mo 4$`d`$ states in this energy region.

The assignment of various spectral features are further confirmed by a comparison of the intensity variations when the photon energy is changed from $`h\nu `$ $`<`$ 100 eV to 1486.6 eV using AlK$`\alpha `$ radiation. In Fig. 3 we show these spectra with the integrated area being normalized for each spectrum, with two representative low photon energy spectra with $`h\nu `$ = 54.6 eV (open circles) and $`h\nu `$ = 86.4 eV (open squares), while the x-ray photoelectron spectrum (XPS, shown with open triangles in Fig. 3) with $`h\nu `$ = 1486.6 eV. It is evident from this comparison that there are substantial changes in the spectral intensities of various features, as the photon energy is changed by $``$ 1400 eV, thereby changing the relative cross-sections of various states drastically. The most obvious spectral changes are the very strong increases in the intensities of features near EF and in the 7-9 eV range, substantial increase in the intensity of the feature at 1.8 eV and a decrease in the energy region 4-6 eV binding energy in the high photon energy XP spectrum. Noting that the Mo 4$`d`$ cross-section increases strongly and that of Fe 3$`d`$ moderately relative to the oxygen 2$`p`$ states, we may attribute substantial Mo 4$`d`$ contributions in the near EF and 7-9 eV binding energy region, primary Fe 3$`d`$ contribution in the 1.8 eV region and a dominance of O 2$`p`$ derived states in the 4-6 eV binding energy. These conclusions are in agreement with the interpretation based on the spectral changes within the low photon energy ($`h\nu `$ $`<`$ 200 eV) range (see Figs. 2(a) and (b)). Fig. 2 also established substantial localized Fe 3$`d`$ contributions in the neighborhood of the 8.3 eV feature, in addition to the Mo 4$`d`$ contribution. These assignments of the spectral features at once provide us with an understanding of both the partial success and the specific failure of the mean-field calculated results (HF-DOS) as well as LDA+U results reported earlier  in the literature, to explain the experimental spectral features. Compared to GGA, the incorporation of Coulomb interaction within the Fe 3$`d`$ manifold in terms of these mean-field approaches improves the agreement between the experiment and the calculation for the 1.8 and 8.3 eV features that have significant Fe 3$`d`$ character (see Fig. 1); this feature also attains substantial Mo 4$`d`$ characters due to large hopping interactions present in the system. However, at this level of mean-field approximation, there is still a rather pronounced disagreement between the experiment and calculation in the 4-6 eV energy range;  our experimental results have already established that this energy region is dominated by O $`p`$ states. Thus, a single-particle approach for the O $`p`$ states fails to describe the experiment, while the avreage effect of Coulomb interactions within Fe 3$`d`$ states improve the description of the Fe 3$`d`$ related features; this suggests a possible role of Coulomb interactions also within the O $`p`$ states. This is not altogether surprising, since it is known from analysis of the KLL Auger spectra for O in transition metal oxides  that the Coulomb interaction strength, $`U_{pp}`$, in the O $`p`$ states is comparable to the Coulomb interaction strength in the transition metal 3$`d`$ states. While $`U_{pp}`$ is not often explicitly used in various approaches to describe the electronic structure of transition metal oxides, same as in the LDA+U method applied so far to Sr2FeMoO6 or the preceeding mean-field results presented in the manuscript, there are some applications of models including the $`U_{pp}`$ term , most notably in the case of copper compounds with substantial $`p`$ admixed character of the charge carriers. Since the charge carriers in Sr2FeMoO6 have significantly admixed Mo 4$`d`$-O 2$`p`$ character , it is reasonable to anticipate a possible role of $`U_{pp}`$ also in the present system, thus suggesting an explanation for the disagreement between the experimental spectra and the calculated results with $`U_{dd}`$ alone. Therefore, we included the $`U_{pp}`$ term in the multiband Hubbard model and calculated the total and partial densities of states within the mean-field approximation; we found that the inclusion of $`U_{pp}`$ = 3.8 eV improves the agreement between the experiment and the calculation substantially, primarily via a renormalization of oxygen related site energies. The agreement improves farther with an enhancement of the splitting of O $`p`$ up- and down-spin bare energies from the value (0.4 eV) estimated from ab initio band structure results, to 0.7 eV; the resulting Fe $`d`$, Mo $`d`$ and O $`p`$ partial DOS are compared with experimental spectra in Fig. 3. Specifically, the dip in the HF-DOS with $`U_{pp}`$ = 0 as well as in GGA results (Fig. 1) in the energy range of 4-6 eV is replaced with a large intensity feature in presence of $`U_{pp}`$, in agreement with the experimentally observed broad maximum. We also find that the partial contributions from Mo 4$`d`$, Fe 3$`d`$ and O 2$`p`$ states in different energy regions are in good agreement with experimental results. Thus, the calculation shows dominant Mo 4$`d`$ contribution due to $`t_{2g}`$ states at $`E_F`$, with significant contributions of $`t_{2g}`$ and $`e_g`$ symmetries also in the energy range of 7-9 eV, dominant Fe 3$`d`$ contributions at 1.5 and 3 eV binding energy regions, arising from $`e_g`$ and $`t_{2g}`$ states respectively, with a significant intensity contribution in the satellite region of 8-9 eV and a dominant O 2$`p`$ nonbonding contribution in the 4-6 eV binding energy.

Spin-polarized DOS along with the total DOS, shown in Fig. 4, establish some interesting features in the electronic structure of this compound. First, the spectral feature appearing near 8.3 eV, outside the GGA calculated bandwidth and associated with correlation-driven satellite, is entirely contributed by the up-spin channel, suggesting a more pronounced localized nature of electron states in this spin channel compared to the down spin channel. This result can be understood in terms of the half-metallic nature of the ground state with the up-spin channel indeed localized with a gap at $`E_F`$, while the down-spin channel is delocalized with a finite DOS at $`E_F`$. However, a gap at the $`E_F`$ does not necessarily imply a localization of the wave-function in space; we confirm the localized-in-space nature of the up-spin channel in contrast to the down-spin channel by plotting the charge densities of the $`e_g`$ and $`t_{2g}`$ spin states close to $`E_F`$ in the basal plane of the unit cell in Fig. 5. This figure shows the extended nature of the down spin charge density with substantial weights on all atomic sites, namely Fe, O and Mo. In clear contrast, the up-spin charge density has virtually no contribution on Mo; consequently, the charge density plot suggests almost completely localized charge densities decoupled and isolated on individual FeO6 octahedra.

We have already pointed out that an analysis of the ab-initio band structure results suggests a splitting of about 0.4 eV between the oxygen up- and down-spin $`p`$ states and that a proper simulation of the experimental spectra required a further enhancement of this splitting to 0.7 eV within the HF approach. This is an unusually large spin-splitting compared to oxygen $`p`$ states in other ferromagnetic oxide systems. One explanation for the spin splitting could be based on the atomic exchange strength. However, since the exchange splitting is given by the product of the moment and the exchange strength and the typical oxygen moment in this compound is less than 0.1 $`\mu _B`$ , this explanation would require an unphysically large exchange strength at the oxygen sites. We have traced back the non-magnetic origin of this spin-splitting to the differing spatial extents of the corresponding oxygen charge densities in up and down spin channels. The overall behavior of the charge densities in the two spin channels are grossly similar and a direct comparison of the two will not show any remarkable difference. In order to enhance the contrasting behaviors between the two spin-channels, we construct the charge density difference between the two spin channels. In Fig. 6 we have plotted this difference in up- and down-spin charge densities in the basal plane. We show the positive part of the $`(updown)`$ charge density in the left panel and the complementary $`(downup)`$ one in the right panel. Thus, the left panel shows the space region with a higher charge density in the up-spin channel compared to the down-spin channel and the right one shows the excess charge density in the down-spin channel compared to the up-spin channel. It is clear from the figure that the down-spin channel of oxygen $`p`$ states are spatially more extended than the up-spin states, with the spatial extension more dominantly in the direction of perpendicular bisector of the Fe-O-Mo bond through the O site. We have analyzed several states within the oxygen non-bonding region; all of them show the common features of more extended down-spin channel with relatively more charge density away from the nearby positively charged metal centers compared to the up-spin channel. This shows that the crystal-field stabilization of the electron states due to the positively charged neighboring sites will be relatively more for the up-spin oxygen states compared to the down-spin ones, as indeed found here, giving rise to the apparent exchange splitting of the oxygen $`p`$ states.

In conclusion, we have found that the electronic structure of Sr2FeMoO6 is dominated by strong electron correlation effects, manifesting more prominently within the Fe 3$`d`$ up-spin states, as a consequence of its half-metallic ground state. Detailed analysis of the spectral features using photon energy variation allows us to establish two interesting effects related to the oxygen $`p`$ states. First, it appears necessary to explicitly account for the Coulomb interactions within the oxygen $`p`$ states and second, there is a splitting between the oxygen $`p`$ up- and down-spin channels arising from different spatial extensions of the corresponding wave-functions. In this context, we note that the ferromagnetic stability in Sr2FeMoO6 is sensitively dependent on the properties of oxygen states and the first principle band structure calculations tend to underestimate the ferromagnetic stability.  Our findings here are likely to provide an understanding for the stability of ferromagnetism in this system, as our preliminary calculation of the spin wave dispersions indicates that the ferromagnetism is indeed favored by both finite $`U_{pp}`$ and the “spin”-splitting of the oxygen $`p`$ states.

Acknowledgement : We thank O. K. Anderson and O. Jepsen for providing the LMTO code. This project is supported by the Department of Science and Technology, Government of India.