text:
# Untitled Document

The Topological Origin of Black Hole Entropy

Dept. of Physics

Beijing Normal University

Beijing 100875, China

Zhong Chao Wu

(Gravity Essay)

Abstract

In gravitational thermodynamics, the origin of a black hole’s entropy is the topology of its instanton or constrained instanton. We prove that the entropy of an arbitrary nonrotating black hole is one quarter the sum of the products of the Euler characteristics of its horizons with their respective areas. The Gauss-Bonnet-like form of the action is not only crucial for the evaluation, but also for the existence of the entropy. This result covers all previous results on the entropy of a nonrotating black hole with a regular instanton. The argument can be readily extended into the lower or higher dimensional model. The problem of quantum creation of such a black hole is completely resolved.

PACS number(s): 98.80.Hw, 98.80.Bp, 04.60.Kz, 04.70.Dy

Keywords: gravitational thermodynamics, quantum cosmology, constrained gravitational instanton, black hole creation

e-mail: wu@axp3g9.icra.it

Hawking radiation was the most important discovery in gravitational physics in the second half of the last century. From this scenario it is derived that the entropy of a Schwarzschild black hole is a quarter of its horizon area . The entropy is interpreted as the measure of our ignorance about the information beyond the horizon. In Euclidean quantum gravity, it was shown that the origin of the entropy is due to the nontrivial topology of its Euclidean spacetime section, i.e., the instanton .

For many cases the concrete relation of the entropy and topology of a black hole has been worked out . It is widely believed that the entropy of a nonrotating black hole is one quarter the sum of the products of the Euler characteristics of its horizons involved and their respective areas. However, this has only been proven for the case of a black hole with a regular instanton. For all known nonrotating black hole cases, a regular instanton can be obtained in three cases: (i) when only one horizon is involved, or (ii) two horizons with the same surface gravities $`\kappa `$ are involved, or (iii) at least one of the two horizons involved recedes into an internal infinity. Only very recently, this relation for an arbitrary nonrotating black hole without a regular instanton has been proven . This allows black holes with regular instantons to be a special case of our considerations here.

In gravitational thermodynamics , it is known that for a generic black hole, one can still define a temperature as $`\kappa /2\pi `$ associated with each horizon. However, there does not exist a thermal equilibrium state with an uniform temperature. Or equivalently, there does not exist a regular instanton. This implies that neither canonical nor grandcanonical ensemble can apply here. Fortunately, one can circumvent this obstacle using microcanonical ensemble. In contrast, for microcanonical ensemble, the temperature is not defined, but all conserved quantities are given. Since the probability of each state under the conserved quantity constraints are equal, the entropy is simply the logarithm of the number of these states.

The partition function for microcanonical ensemble in gravitational thermodynamics is the path integral

$$Z=d[g_{\mu \nu }]d[\varphi ]\mathrm{exp}I$$
(1)
where the path integral is over all spacetime metrics $`g_{\mu \nu }`$ and matter configurations $`\varphi `$ under the conditions for the ‘equator’ imposed by the restriction of the microcanonical ensemble, and $`I`$ is the Euclidean action.

The Euclidean action in the Einstein theory is

$$I=\frac{1}{16\pi }_M\sqrt{g}(R2\mathrm{\Lambda }+L_m)d^4x\frac{1}{8\pi }_M\sqrt{g}Kd^3x,$$
(2)
where $`R`$ is the scalar curvature of the spacetime manifold $`M`$, $`K`$ is the expansion rate of the boundary $`M`$, $`g`$ is the determinant of the metric for the 4-metric or its lower dimensional version, $`\mathrm{\Lambda }`$ is the cosmological constant, and $`L_m`$ is the Lagrangian of the matter content.

The main contribution to the path integral is from a stationary action orbit. The partition function is approximated by the exponential of the negative of the action of the orbit. This is called the $`WKB`$ approximation, which we shall use in this paper.

If the action of the orbit is stationary with respect to all variations, then one obtains an instanton. It is determined essentially by the topological properties of the manifold. The metric is regular without any singularity.

On the other hand, the dominating contribution to the partition function for the microcanonical ensemble is called constrained instanton . It is a manifold for which the action is stationary with respect only to the variations under the restrictions due to the ensemble, instead of with respect to all variations under no restriction.

For the quantum creation scenario in the no-boundary universe , the relative probability takes the same form (1) and the path integral is over all 4-metrics with the given 3-metric and matter content at the equator. These constraints can be characterized by a few parameters, like mass $`m`$, charge $`Q`$ and angular momentum $`J`$ for the black hole case. These conditions are exactly the same as for the microcanonical ensemble in gravitational thermodynamics. Therefore, the constrained instanton is also the creation seed in the no-boundary universe. The exponential of the negative of the instanton action is the relative creation probability of the universe.

Since the entropy is the logarithm of the partition function in microcanonical ensemble, then at the same level, the entropy is the negative of the action.

For all cases of black holes considered, the spacetime has a $`U(1)`$ isometry. The group parameter is identified as the Killing time coordinate. The Euclidean nonrotating black hole metric takes the form

$$ds^2=\mathrm{\Delta }(r)d\tau ^2+\mathrm{\Delta }^1(r)dr^2+r^2d\mathrm{\Omega }_2^2,$$
(3)
where $`\tau =it`$, and the 2-metric $`d\mathrm{\Omega }_2^2`$ is a compact manifold which does not depend on coordinates $`\tau `$ and $`r`$. For an ordinary black hole the 2-metric is a sphere, and for a topological black hole it is a compactified plane or hyperboloid. The zeros of the rational function $`\mathrm{\Delta }(r)`$ are the horizons.

One can construct a compact constrained instanton (to be justified below) using a sector between two horizons denoted by two zeros $`r_l,r_k`$ in the identified manifold. The surface gravity $`\kappa _i`$ of the horizon $`r_i`$ is $`d\mathrm{\Delta }(r)/2dr|_{r=r_i}`$. If the zero is of multiplicity 1, then one obtains a nonzero $`\kappa _i`$. On the two dimensional space $`(\tau ,r)`$ the conical singularity at the horizon can be regularized by choosing $`\beta =2\pi \kappa _i^1`$. For two horizons with same nonzero surface gravities, one can use the same $`\beta `$ to obtain a compact regular instanton. If these two surface gravities are different, then the constrained instanton has at least one conical singularity at the horizons, since no value of parameter $`\beta `$ can regularize both of the horizons simultaneously. The de Sitter model is an exception, since $`r`$ is identified with $`r`$, only one horizon, i.e. the cosmological horizon is needed for the construction of the instanton.

If one of the two zeros is of multiplicity larger than 1, then its surface gravity $`\kappa _i=0`$, and the associated horizon recedes to an internal infinity. Then it is always possible to regularize the other horizon by choosing a right value $`\beta `$ to obtain a regular instanton. The most familiar case is the extreme Reissner-Nordstr$`\ddot{𝐨}`$m black hole in the nonvacuum model .

Now, let us calculate the action of the pasted manifold. We use $`M_l`$ to denote the small neighbourhood of horizon $`r_l`$ with the boundary of a constant coordinate $`r`$. The Euler number $`\chi (l)`$ for the 2-dimensional $`(\tau ,r)`$ section of neighbourhood with zero (nonzero) surface gravity is 0 (1). We use $`M^{}`$ to denote $`M`$ minus $`M_l`$ and $`M_k`$. For the form of action (2), the total action is the sum of those from the three submanifolds.

First of all, let us consider the vacuum model with a cosmological constant. The total action is

$$I=I_l+I_k+_M^{}(\pi ^{ij}\dot{h}_{ij}NH_0N_iH^i)d^3x𝑑\tau ,$$
(4)
where the actions $`I_l`$ and $`I_k`$ are the actions for $`M_l`$ and $`M_k`$. The action of $`M^{}`$ has been recast into the canonical form. $`N`$ and $`N_i`$ are the lapse function and shift vector, $`h_{ij}`$ and $`\pi ^{ij}`$ are the 3-metric and the conjugate momenta respectively, $`H_0`$ and $`H^i`$ are the Einstein and momentum constraints, and the dot denotes the time derivative. The manifold satisfies the Einstein equation, and all time derivatives vanish due to the $`U(1)`$ isometry, therefore the integral over $`M^{}`$ is equal zero.

Now the action $`I_l`$ or $`I_k`$ can be written

$$I_i=\frac{1}{16\pi }_{M_i}\sqrt{g}(R2\mathrm{\Lambda })d^4x\frac{1}{8\pi }_{M_i}\sqrt{g}Kd^3x,(i=l,k).$$
(5)
If there is a conical singularity at the horizon, its contribution to the action can be considered as the degenerate version of the second term of the action, in addition to that from the boundary of $`M_i`$. The conical singularity contribution is termed as a deficit “angle” due to its emergence. If the horizon recedes into an internal infinity, then this is no longer of concern.

One can apply the Gauss-Bonnet theorem to the 2-dimensional $`(\tau ,r)`$ section of $`M_i`$,

$$\frac{1}{4\pi }_{\widehat{M}_i}\sqrt{g}^2Rd^2x+\frac{1}{2\pi }_{\widehat{M}_i}\sqrt{g}^1Kd^1x+\frac{\delta _i}{2\pi }=\chi (i),$$
(6)
where $`\widehat{M}_i`$ is the projection of $`M_i`$ onto the 2-dimensional $`(\tau ,r)`$ section, $`{}_{}{}^{2}R`$ is the scalar curvature on it, $`{}_{}{}^{1}K`$ is the corresponding expansion rate, $`\delta _i`$ is the total deficit angle, and $`\chi (i)`$ is the Euler characteristic of $`\widehat{M}_i`$. Since the expansion rate of the subspace $`r^2d\mathrm{\Omega }_2^2`$ goes to zero at the horizon, $`K`$ and $`{}_{}{}^{1}K`$ are equal. Comparing eqs (5) and (6), one can see that as the circumference of the boundary tends to infinitesimal, the action (5) becomes $`\chi (i)A_i/4`$, where $`A_i`$ is the surface area of the horizon. It is noted that both the volume integral of (5) and the first term of the left hand side of (6) vanish as the boundary approaches the horizon. The same result is obtained regardless of the existence of the conical singularity at the horizon or not, i.e., it is independent of the value $`\beta `$. Here the conical singularity contribution is represented by $`\delta _i/2\pi `$.

From (4)-(6), we learn that the action is independent of the parameter $`\beta `$. Since the manifold satisfies the Einstein equation everywhere with probable exception at the horizons. The conical singularities there are parametrized by $`\beta `$, the only degree of freedom left. Therefore, the manifold is qualified as a constrained instanton. The entropy, or the negative of the total action of the constrained instanton is

$$S=I=\frac{1}{4}(\chi (l)A_l+\chi (k)A_k).$$
(7)
This is a quite universal formula.

If one use $`r_0`$ to denote the maximum zero of $`\mathrm{\Delta }(r)`$ and the metric of the sector $`r>r_0`$ is Euclidean, then this sector can also be used for constructing an open instanton.

The action of the open instanton is divergent. One can modify the action form (4) to regularize it, by setting $`I_l=I_0`$, then dropping $`M_k`$ term and finally letting $`M^{}`$ be the sector $`M`$ minus $`M_0`$. The form of this action is derived from the requirement that, as the mass $`m`$ (and the angular momentum $`J`$ for a rotating model not expressed by (3)) is held fixed at infinity with an appropriate asymptotic falloff for the field, the Einstein and field equations must be derived from the action . These boundary conditions correspond to the microcanonical ensemble. Follow the same argument, one can derive $`S=\chi (0)A_0/4`$.

One can include the Maxwell or gauge field into the model. For the magnetic black hole, the Maxwell action is compatible with the requirement of the microcanonical ensemble. On the other hand, for the electric black hole, it is not. One has to appeal to a Legendre transformation here. In both cases, formula (7) remains valid .

It is noted that the action of the black hole metric must be a linear function of the parameter $`\beta `$ due to the $`U(1)`$ isometry. The independence of the action from the parameter is at the edge of the knife. Indeed, the Gauss-Bonnet-like form is crucial not only for the evaluation, but also for the existence of the entropy! The Equivalence Principle is not sufficient to restrict the action to take the form (2) in the 4-dimensional spacetime. The deep reason behind this remains a mystery in Nature.

The method of the dimensional continuation of the Gauss-Bonnet theorem has been used to study the entropy of a black hole with a regular instanton . The formula presented in  does not apply to a topological black hole . In contrast, formula (7) is true for all nonrotating black holes.

Our analysis can be generalized into the Lovelock theory of gravitation . One can study $`n`$dimensional black holes, which are described by (3) with $`d\mathrm{\Omega }_2^2`$ replaced by $`d\mathrm{\Omega }_{n2}^2`$. The entropy of a black hole or the negative of the action is

$$S=I=\frac{1}{4}(\chi (l)A_lf_l+\chi (k)A_kf_k),$$
(8)
where $`f_i`$ is a numerical factor determined by the metric $`r^2d\mathrm{\Omega }_{n2}^2`$ of the horizon.

The discussion can be extended straightforward to the nonvacuum model. The formula for the entropy of a black hole with a regular instanton has been previously obtained . The constrained instanton can be used as a seed for quantum creation of a black hole in Einstein gravity and its higher or even lower dimensional black hole . However, one has to use the instanton with the largest action . After we have obtained formulas (7) and (8), we no longer need to check stationary property of the constrained instanton case by case for an arbitrary black hole described by (3) and its lower or higher dimensional version. The entropy and quantum creation of an arbitrary nonrotating black hole are thus completely resolved.

References:

1. S.W. Hawking, Commun. Math. Phys. 43, 199 (1975).

2. S.W. Hawking, in General Relativity: An Einstein Centenary Survey, eds. S.W. Hawking and W. Israel, (Cambridge University Press, 1979).

3. G.W. Gibbons and R.E. Kallosh, Phys. Rev. D51, 2839 (1995).

4. S.W. Hawking, G.T. Horowitz and S.F. Ross, Phys. Rev. D51, 4302 (1995).

5. Z.C. Wu, gr-qc/9911078, (in press).

6. G.W. Gibbons and S.W. Hawking, Phys. Rev. D15, 2752 (1977).

7. G.W. Gibbons and S.W. Hawking, Phys. Rev. D15, 2738 (1977).

8. Z.C. Wu, Int. J. Mod. Phys. D6, 199 (1997), gr-qc/9801020.

9. Z.C. Wu, Phys. Lett. B445, 274 (1999), gr-qc/9810012.

10. J.B. Hartle and S.W. Hawking, Phys. Rev. D28, 2960 (1983).

11. C. Teitelboim, Phys. Rev. D51, 4315 (1995).

12. S.W. Hawking and G.T. Horowitz, Class. Quant. Grav. 13, 1487 (1996).

13. T. Regge and C. Teitelboim, Ann. Phys. (N.Y.) 88, 286 (1974).

14. M. Ba$`\stackrel{~}{n}`$ados, C. Teitelboim and J. Zanelli, Phys. Rev. Lett. 72, 957 (1994).

15. S. Liberati and G. Pollifrone, Phys. Rev. D56, 6458 (1997).

16. D.R. Brill, J. Louko and P. Peldan, Phys. Rev. D56, 3600 (1997).

17. Z.C. Wu, gr-qc/9907064.

18. D. Lovelock, J. Math. Phys. 12, 498 (1971).

19. C. Teitelboim and J. Zanelli, Class. Quant. Grav. 4, L125 (1987).

20. T. Jacobson and R. Myers, Phys. Rev. Lett. 70, 3684 (1993).

21. R.C. Myers and J.Z. Simon, Phys. Rev. D38, 2434 (1988).

22. Z.C. Wu, gr-qc/9907065.


format: text

text:
# Lattice Points in Lattice Polytopes

## 1 Introduction

A lattice polytope in $`^d`$ is a convex polytope whose vertices are lattice points, that is, points in $`^d`$. For an integer $`l1`$, let $`I_l(P)=\text{int}(P)l^d`$ be the set of interior points of $`P`$ whose coordinates are integers divisible by $`l`$.

Of course, some points of $`I_l(P)`$ can lie ‘close’ to $`P`$, the boundary of $`P`$. However, our Theorem 4 shows that, provided $`I_l(P)\mathrm{}`$, there is $`𝐰I_l(P)`$ with

$$\text{ca}(𝐰,P)8d(8l+7)^{2^{2d+1}},$$
(1)
where $`\text{ca}(𝐰,P)`$ is the coefficient of asymmetry of $`P`$ about $`𝐰`$:

$$\text{ca}(𝐰,P)=\underset{|𝐲|=1}{\mathrm{max}}\frac{\mathrm{max}\{\lambda 𝐰+\lambda 𝐲P\}}{\mathrm{max}\{\lambda 𝐰\lambda 𝐲P\}}.$$

Although the function in the right-hand side of (1) is enormous, the main point is that it depends only on $`d`$ and $`l`$.

We prove an inequality of this type for the case of a simplex $`S`$ first. Namely, Theorem 2 implies that, for some $`𝐰I_l(S)`$,

$$\text{ca}(𝐰,S)8(8l+7)^{2^{d+1}}.$$
(2)
Here the claim essentially concerns the barycentric coordinates $`(\alpha _0,\mathrm{},\alpha _d)`$ of $`𝐰`$ inside $`S`$ because of the easy relation

$$\text{ca}(𝐰,S)=\underset{0id}{\mathrm{max}}\frac{1\alpha _i}{\alpha _i}=\frac{1}{m_S(𝐰)}1,$$
(3)
where $`m_S(𝐰):=\mathrm{min}_{0id}\alpha _i`$ is the smallest barycentric coordinate of $`𝐰S`$. Define

$$\beta (d,l):=\underset{S}{inf}\mathrm{max}\{m_S(𝐰)𝐰I_l(S)\},$$
(4)
where the infimum is taken over all lattice simplices $`S`$ with $`I_l(S)\mathrm{}`$. (For example, it is easy to see that $`\beta (1,l)=\frac{1}{l+1}`$.) Thus we have to prove a positive lower bound on $`\beta (d,l)`$. The gist the proof is that if we have $`𝐰I_l(S)`$ with $`m_S(𝐰)`$ being ‘small’, then using one approximation lemma of Lagarias and Ziegler  we can ‘jump’ to another vertex $`𝐰^{}I_l(S)`$ with $`m_S(𝐰^{})>m_S(𝐰)`$, see Theorem 2.

In fact, one result of Lawrence \[4, Theorem 3\] implies that $`\beta (d,1)>0`$ but does not give any explicit bound, see Section 7 here.

It would be interesting to know how far our bounds (1) and (2) are from the best possible values. The best values that we know arise from the following family of lattice simplices.

Define inductively the sequence $`t_{d,l}`$ by $`t_{1,l}=l+1`$ and $`t_{d+1,l}=t_{d,l}^2t_{d,l}+1`$. (This sequence appears in .) Consider the simplex

$$B_{d,l}:=\text{conv}\{t_{1,l}𝐞_1,\mathrm{},t_{d1,l}𝐞_{d1},t_{d,l}𝐞_d,𝐞_d\}^d,$$
(5)
where $`(𝐞_1,\mathrm{},𝐞_d)`$ is the standard basis. It is not hard to see that $`I_l(B_{d,l})=\{l\mathrm{𝟏},l(\mathrm{𝟏}𝐞_d)\}`$, cf. \[3, Proposition 2.6\]. We have $`m_{B_{d,l}}(l\mathrm{𝟏})=m_{B_{d,l}}(l(\mathrm{𝟏}𝐞_d))=\frac{l}{t_{d,l}^21}`$: the vertex $`l\mathrm{𝟏}=(l,\mathrm{},l)`$, for example, has barycentric coordinates

$$(\frac{l}{t_{1,l}},\mathrm{},\frac{l}{t_{d1,l}},\frac{l}{t_{d,l}1}\times \frac{t_{d,l}}{t_{d,l}+1},\frac{l}{t_{d,l}1}\times \frac{1}{t_{d,l}+1}).$$

One can show that $`t_{d,l}(l+1)^{2^{d2}}+1`$ for $`d2`$ by considering $`u_{d,l}=t_{d,l}1`$; hence

$$\beta (d,l)\frac{l}{t_{d,l}^21}l(l+1)^{2^{d1}},d2.$$
(6)
Thus (2) establishes the correct type of dependence on $`d`$ and $`l`$, although the gap between the bounds is huge. Perhaps, $`B_{d,l}`$ gives the actual value of the function $`\beta (d,l)`$ as well as the sharp bound for (1).

To extend Theorem 2 to a general lattice polytope $`P^d`$ (Theorem 4), we try to find a lattice polytope $`P^{}P`$ with few vertices such that $`I_l(P^{})\mathrm{}`$ and a homothetic copy of $`P^{}`$ covers $`P`$. The latter condition gives an upper bound on $`\text{ca}(𝐰,P)`$ in terms of $`\text{ca}(𝐰,P^{})`$ for $`𝐰\text{int}(P^{})`$, see Lemma 3, and is satisfied if, for example, $`P^{}S`$, where $`SP`$ is a simplex of the maximum volume. But to get a non-empty $`I_l(P^{})`$ we may have to add as many as $`d`$ extra vertices to $`S`$. It is now possible to define our jumps within $`P^{}`$ to get the required $`𝐰I_l(P^{})`$. However, the bound (1) for $`d`$-polytopes that we obtain is comparable with that for $`2d`$-dimensional simplices; we believe that we lose here too much but we have not found any better argument.

Next, we investigate the following problem. Let $`p(d,k,l)`$ (resp. $`s(d,k,l)`$) be the maximum volume of a lattice polytope (resp. simplex) $`P^d`$ with $`|I_l(P)|=k`$. As for any $`d2`$ there exist lattice simplices with no lattice points in the interior and of arbitrarily large volume, we restrict our consideration to the case $`k1`$.

Trivially, $`p(1,k,l)=s(1,k,l)=(k+1)l`$. A result of Scott  implies that $`p(2,1,1)=s(2,1,1)=\frac{9}{2}`$ and $`p(2,k,1)=s(2,k,1)=2(k+1)`$ for $`k2`$. Hensley \[2, Theorem 3.6\] showed that $`p(d,k,1)`$ exists (i.e., it is finite) for $`k1`$. The method of Hensley was sharpened by Lagarias and Ziegler \[3, Theorem 1\], who showed that

$$p(d,k,l)kl^d(7(kl+1))^{d2^{d+1}},$$
(7)
and also observed that, for any fixed $`(d,k,l)`$, there are finitely many (up to a $`GL_n()`$-equivalence) lattice polytopes $`P^d`$ with $`|I_l(P)|=k1`$.

Lagarias and Ziegler \[3, Theorem 2.5\] proved the following extension of a theorem of Mahler : “A convex body $`K^d`$ with $`k=|I_l(K)|1`$ satisfies

$$\text{vol}(K)\left(l(\text{ca}(𝐰,K)+1)\right)^dk,$$
(8)
for any $`𝐰I_l(K)`$.”

Combining (8) with (1) (or more exactly with (27)), we obtain that

$$p(d,k,l)(8dl)^d(8l+7)^{d2^{2d+1}}k.$$
(9)
A theorem of Blichfeldt  says that $`|P^d|n+n!\text{vol}(P)`$; combined with (9) it gives an upper bound on $`|P^d|`$ in terms of $`|I_l(P)|`$ (if the latter set is non-empty).

An upper bound on $`s(d,k,l)`$ can be obtained by applying (8) to (2). However, we obtain a better bound in Theorem 6 by exploiting the geometry of a simplex, namely we show that

$$s(d,k,l)2^{3d2}l^d(8l+7)^{(d1)2^{d+1}}k/d!.$$
(10)
The best lower bound on $`p(d,k,l)`$ and $`s(d,k,l)`$ that we know (except for $`(d,k,l)=(2,1,1)`$), comes from the consideration of the simplex

$$S_{d,k,l}:=\text{conv}\{\mathrm{𝟎},t_{1,l}𝐞_1,\mathrm{},t_{d1,l}𝐞_{d1},(k+1)(t_{d,l}1)𝐞_d\},$$

which satisfies $`I_l(S_{d,k,l})=\{l\mathrm{𝟏}+il𝐞_d0ik1\}`$, see \[3, Proposition 5.6\]. This demonstrates that

$$s(d,k,l)\text{vol}(S_{d,k,l})>\frac{k+1}{d!l}(l+1)^{2^{d1}},$$

see formula (2.13) in . The family $`(S_{d,k,1})`$ was found by Zaks, Perles and Wills  and its generalization (the addition of parameter $`l`$)—by Lagarias and Ziegler .

Again, we have the correct type of dependence of $`d`$, $`k`$ and $`l`$ but the gap between the known bounds is huge. The ultimate aim would be to find exact values, which is probably not hopeless because the above contructions, believed to be extremal, are rather simple.

## 2 Jumping inside a simplex

We will use the following lemma of Lagarias and Ziegler \[3, Lemma 2.1\].

###### Lemma 1

For a real $`\lambda 1`$ and integer $`n1`$, define

$$\delta (n,\lambda )=\left(7(\lambda +1)\right)^{2^{n+1}}.$$
(11)
Then, for all positive real numbers $`\alpha _1,\mathrm{},\alpha _n`$ satisfying

$$1\delta (n,\lambda )<\underset{i=1}{\overset{n}{}}\alpha _i1,$$

there exist non-negative integers $`P_1,\mathrm{},P_n,Q`$ such that

$`Q`$ $`=`$ $`P_1+\mathrm{}+P_n>0,`$ (12)
$`\alpha _i`$ $`>`$ $`{\displaystyle \frac{\lambda P_i}{\lambda Q+1}}\text{ for }1in,`$ (13)
$`\lambda Q+1`$ $``$ $`\delta (n,\lambda )^1.\text{ }\text{ }`$ (14)
The above lemma is the main ingredient in our ‘jumps.’ Here it is applied with $`\lambda =\frac{8}{7}l`$. There is nothing special about the constant $`\frac{8}{7}`$ except it makes (11) look simpler; any fixed number greater than $`1`$ would do as well.

###### Theorem 2

Let $`l1`$ and let $`S=\text{conv}\{𝐯_0,\mathrm{},𝐯_d\}^n`$ be a lattice simplex. If $`\text{rel-int}(S)l^n`$ is non-empty, then it contains a point $`𝐰`$ with

$$m_S(𝐰)\gamma :=\delta (d,\frac{8}{7}l)/8=(8l+7)^{2^{d+1}}/8.$$
(15)
Proof. We may assume that $`n=d`$ because we can always find a linear transformation preserving the lattice $`^n`$ (and so $`l^n`$ as well) and mapping $`S`$ into $`^d^n`$.

Let $`𝐰=_{i=0}^d\alpha _i𝐯_iI_l(S)`$, $`_{i=0}^d\alpha _i=1`$, be a vertex maximizing $`m_S(𝐰)`$. Suppose that the claim is not true. Assume that $`\alpha _0\mathrm{}\alpha _d`$; then $`m_S(𝐰)=\alpha _0<\gamma `$. Let $`j`$ be the index with $`\alpha _j<8\gamma \alpha _{j+1}`$; note that $`jd1`$ is well-defined as $`\alpha _d\frac{1}{d+1}8\gamma `$.

We have $`_{i=0}^j\alpha _j<8\gamma (j+1)`$ which, as it is easy to see, does not exceed $`\delta (dj,\frac{8}{7}l)`$ for $`j[0,d1]`$. Hence, Lemma 1 is applicable to the $`dj`$ numbers $`\alpha _{j+1},\mathrm{},\alpha _d`$ and yields integers $`P_{j+1},\mathrm{},P_d,Q`$ satisfying (12)(14).

Consider the vertex

$$𝐰^{}=(lQ+1)𝐰\underset{i=j+1}{\overset{d}{}}lP_i𝐯_il^d.$$

We have $`𝐰^{}=_{i=0}^d\alpha _i^{}𝐯_i`$, where, for $`i[0,j]`$, $`\alpha _i^{}:=(lQ+1)\alpha _i>\alpha _0`$ and, for $`i[j+1,d]`$, $`\alpha _i^{}:=(lQ+1)\alpha _ilP_i>\alpha _i/8\alpha _0`$ by (13). As

$$\underset{j=0}{\overset{n}{}}\alpha _j^{}=(lQ+1)\underset{i=0}{\overset{d}{}}\alpha _i\underset{i=j+1}{\overset{d}{}}lP_i=1,$$

the lattice point $`𝐰^{}`$ lies in the interior of $`S`$ and contradicts the choice of $`𝐰`$.

Remark. For $`n=d`$ the inequality (2) claimed in the introduction follows by applying (3) to the vertex $`𝐰I_l(S)`$ given by Theorem 2

## 3 $`\beta (d,l)`$ for small $`d`$ and $`l`$

Let us try to deduce some estimates of $`\beta (d,l)`$ when $`d`$ and $`l`$ are small. We have a general upper bound (6) which in particular says that

$`\beta (2,l)`$ $``$ $`{\displaystyle \frac{l}{t_{2,l}^21}}={\displaystyle \frac{1}{(l+1)(l^2+l+2)}},`$ (16)
$`\beta (3,1)`$ $``$ $`{\displaystyle \frac{1}{48}},\beta (3,2){\displaystyle \frac{1}{924}}.`$ (17)
Here we present some results obtained with the help of computer showing that (16) and (17) are probably sharp.

How can one get a lower bound on, for example, $`\beta (2,1)`$? Our approach is the following.

Given a lattice simplex $`S=\text{conv}\{𝐯_0,𝐯_1,𝐯_2\}^2`$, let $`𝐰`$ be a lattice vertex maximizing $`m_S`$ over $`I_1(S)\mathrm{}`$. Write the barycentric representation $`𝐰=_{i=0}^2\alpha _i𝐯_i`$. We have

$$\alpha _0+\alpha _1+\alpha _2=1.$$
(18)
Without loss of generality we may assume that

$$0\alpha _0\alpha _1\alpha _2.$$
(19)
Consider the vertex $`𝐰^{}=2𝐰𝐯_3`$ which is the jump of $`𝐰`$ corresponding to $`𝐏=(0,0,1)`$. Its barycentric coordinates $`(\alpha _0^{},\alpha _1^{},\alpha _2^{})=(2\alpha _0,2\alpha _1,2\alpha _21)`$ satisfy $`\alpha _1^{}\alpha _0^{}>\alpha _0`$. By the choice of $`𝐰`$, we must have

$$2\alpha _21\alpha _0.$$
(20)
Similarly, the $`(0,1,1)`$-jump $`𝐰^{}=3𝐰𝐯_1𝐯_2`$ satisfies $`\alpha _0^{}=3\alpha _0>\alpha _0`$ and, in the view of $`\alpha _2^{}\alpha _1^{}`$, we obtain

$$3\alpha _11\alpha _0.$$
(21)
Not everything goes so smoothly if we consider e.g. the $`(0,1,2)`$-jump, when we can only deduce that

$$4\alpha _11\alpha _0\text{or}4\alpha _22\alpha _0.$$
(22)
How small can $`\alpha _0`$ be, given only the constraints (18)(22) (which can be realized as a mixed integer program)? Solving this MIP, we obtain that $`\alpha _0\frac{2}{19}`$. Knowing this bound, we can enlarge our arsenal of jumps. For example, the vertex $`𝐰^{}=12𝐰𝐯_14𝐯_26𝐯_2`$ satisfies $`\alpha _0^{}=12\alpha _01>\alpha _0`$; hence

$$12\alpha _14\alpha _0\text{or}12\alpha _26\alpha _0.$$
(23)
The addition of (23) to the system (18)(22), improves our lower bound to $`\alpha _0\frac{2}{17}`$.

In this manner we can repeatedly add new constraints to our MIP as long as this improves the lower bound on $`\alpha _0`$. This was realized as a program in C which can be linked with either CPLEX (commercial) or lp\_solve (public domain) MIP solver. The source code is freely available  and the reader is welcome to experiment with it.

When one runs the program, it seems that the obtained lower bound $`f`$ on $`a_0`$ approaches some limit without attaining it. And, of course, the more iterations we perform, the larger coefficients in the added inequalities are and the problem becomes more complex. Hence, the question to what extend we can trust the output should be considered. The CPLEX has the mechanism to set up various tolerances which specify how far CPLEX allows variables to violate the bounds and to be still considered feasible during perturbations. The manual does not specify the guaranteed accuracy. As the coefficient at $`\alpha _0`$ is always $`1`$, we assume that the error does not exceed $`\mathrm{\Delta }`$, the largest sum of absolute values of coefficients in one inequality multiplied by the tolerance, which we set to $`10^9`$, the smallest value allowed by CPLEX’s manual.

In particular, we allow a $`(P_0,P_1,P_2)`$-jump when we can guarantee that $`(1+P_0+P_1+P_2)\alpha _0P_0\alpha _0`$, namely, when $`(f\mathrm{\Delta })(P_0+P_1+P_2)P_0`$.

We ran the program, with CPLEX 6.6, for various $`d`$ and $`l`$; Table 1 records the obtained lower bounds up to $`10^6`$. Unfortunately, we had no success when $`d4`$ or when $`d=3`$ and $`l3`$ or when $`d=2`$ and $`l28`$: the obtained lower bound was still zero when the MIP became too large to solve.

## 4 Extending results to lattice polytopes

First we have to express analytically the intuitively obvious fact that if two polytopes cover each other (up to a small homothety) then their coefficients of asymmetry cannot be far apart.

###### Lemma 3

Let $`P^{}P`$ be two polytopes such that $`P`$ can be covered by a translate of $`\lambda P^{}`$. Then, for any $`𝐰\text{int}(P^{})`$,

$$\text{ca}(𝐰,P)|\lambda |\text{ca}(𝐰,P^{})+|\lambda |1.$$
(24)
Proof. Assume that $`|\lambda |>1`$, for otherwise $`P^{}=P`$ and we are home. Also, the case of $`1`$-dimensional polytopes is trivial.

Let $`𝐰_1,𝐰_2P`$ be two points with $`𝐰[𝐰_1,𝐰_2]`$ and $`\text{ca}(𝐰,P)=|𝐰_1𝐰|:|𝐰𝐰_2|`$. Let $`P^{}[𝐰_i,𝐰]=\{𝐰_i^{}\}`$, $`i=1,2`$, where $`[𝐱,𝐲]`$ denotes the straight line segment between $`𝐱`$ and $`𝐲`$. Clearly,

$$|𝐰_1𝐰_2|=(\text{ca}(𝐰,P)+1)|𝐰𝐰_2|(\text{ca}(𝐰,P)+1)|𝐰𝐰_2^{}|.$$
(25)
As $`\lambda P^{}`$ covers $`\{𝐰_1,𝐰_2\}`$, there are $`𝐮_1,𝐮_2P^{}`$ with

$$𝐰_1𝐰_2=|\lambda |(𝐮_1𝐮_2).$$
(26)
We can assume that $`𝐮_1,𝐮_2P^{}`$. If $`𝐮_1=𝐰_1^{}`$ and $`𝐮_2=𝐰_2^{}`$, then we let $`𝐯=𝐮_2`$. Otherwise let $`𝐯`$ be the (well-defined) point of intersection of $`L(𝐮_1,𝐰)`$ and $`L(𝐮_2,𝐰_2^{})`$, where $`L(𝐱,𝐲)`$ denotes the line though the points $`𝐱`$ and $`𝐲`$. As $`𝐯L(𝐮_2,𝐰_2^{})`$ lies outside of $`\text{int}(P^{})`$, we have

$$\text{ca}(𝐰,P^{})\frac{|𝐮_1𝐰|}{|𝐰𝐯|}=\frac{|𝐮_1𝐮_2||𝐰𝐰_2^{}|}{|𝐰𝐰_2^{}|},$$

which implies the required by (25) and (26).

Remark. Note that the bound in (24) is sharp, as is demonstrated e.g. by $`C_d|\lambda |C_d`$ and $`𝐰=c\mathrm{𝟏}`$ with $`0<c\frac{1}{2}`$, where $`C_d^d`$ is the $`0/1`$-cube.

Now we are ready to prove our result on lattice polytopes.

###### Theorem 4

Let $`l1`$ be an integer and let $`P^d`$ be a lattice polytope with $`I_l(P)\mathrm{}`$. Then there is $`𝐰I_l(P)`$ with

$$\text{ca}(𝐰,P)\frac{8d}{\delta (2d,\frac{8}{7}l)}1=8d(8l+7)^{2^{2d+1}}1.$$
(27)
Proof. Let $`S=\text{conv}\{𝐯_0,\mathrm{}𝐯_d\}P`$ be a simplex of the maximum volume; we may assume that each $`𝐯_i`$ is a vertex of $`P`$.

Choose $`𝐮I_l(P)`$. Let $`𝐮_1\text{int}(S)`$ be any vertex and let $`𝐮_2`$ be the point of intersection of the ray $`\{𝐮_1+\lambda (𝐮𝐮_1)\lambda 0\}`$ with the boundary of $`P`$. The vertex $`𝐮_2`$ lies in the interior of some face which is spanned by at most $`d`$ vertices of $`P`$. Hence $`𝐮\text{rel-int}([𝐮_1,𝐮_2])`$ can be represented as a positive convex combination of $`n2d+1`$ vertices of $`P`$ including all vertices of $`S`$, let us say $`𝐮=_{i=0}^n\alpha _i^{}𝐯_i`$ with $`_{i=0}^n\alpha _i^{}=1`$ and each $`\alpha _i^{}>0`$.

Choose a vertex $`𝐰I_l(P^{})`$ and a representation $`𝐰=_{i=0}^n\alpha _i𝐯_i`$ with $`_{i=0}^n\alpha _i=1`$ maximizing $`\mathrm{min}_{0in}\alpha _i`$. Denote this maximum by $`m(𝐰)>0`$. The argument of Theorem 2 shows that $`m(𝐰)\delta (n1,\frac{8}{7}l)/8\delta (2d,\frac{8}{7}l)/8`$.

The polytope $`P^{}`$ can be represented as a projection of an $`n`$-simplex $`S_n`$ such that $`𝐰`$ is the image of $`𝐯\text{int}(S_n)`$ with $`m_{S_n}(𝐯)=m(𝐰)`$. Now, it is easy to see that a linear mapping cannot increase the coefficient of asymmetry; hence

$$\text{ca}(𝐰,P^{})\text{ca}(𝐯,S_n)=\frac{1m(𝐰)}{m(𝐰)}.$$

It is known that $`P(d)S+(d+1)𝐬`$, where $`𝐬`$ is the centroid of $`S`$, see e.g. \[3, Theorem 3\]. By Lemma 3, we obtain

$$\text{ca}(𝐰,P)d\text{ca}(𝐰,P^{})+d1d\frac{1m(𝐰)}{m(𝐰)}+d1=\frac{d}{m(𝐰)}1,$$

which give the required by (11).

Remark. The bound (27) is much worse than (2); the reason is that we may have to approximate $`2d`$-tuples of numbers in Lemma 1. Unfortunately, we cannot guarantee that $`P^{}`$ has much fewer than $`2d+1`$ vertices, as e.g. $`P=\text{conv}\{\pm 𝐞_1,\mathrm{},\pm 𝐞_{d1},(k+1)l𝐞_d\}`$ demonstrates. Perhaps, one can show that any such example cannot be extremal for our problem and thus improve on (27).

Remark. It should be possible to generalize Theorem 2 and 4 by proving the existence of a number $`b=b(d,l,m)>0`$ such that any lattice polytope $`P`$ contains $`m`$ distinct points in $`I_l(P)`$ (provided $`|I_l(P)|m`$) with coefficient of asymmetry of each being at least $`b`$. The idea of the proof is the following. If $`P`$ is a simplex, take distinct $`𝐰_1,\mathrm{},𝐰_mI_l(P)`$ with with largest $`m_P`$’s. Now each jump of $`𝐰_i`$ either does not increase $`m_P(𝐰_i)`$ or maps $`𝐰_i`$ into some other $`𝐰_j`$. We are done if we can show that if $`m_P`$ is very small then there are at least $`m`$ distinct jumps increasing it. The latter claim would be achieved by rewriting the proof of Lemma 1, so that in the conclusion we have at least $`m`$ suitable $`(n+1)`$-tuples of integers. To extend the claim to general lattice polytopes, observe that $`m`$ vertices in $`I_l(P)`$ can be represented each as a positive combination of $`d+1`$ vertices of a max-volume symplex and at most $`md`$ other vertices of $`P`$ and follow the argument of Theorem 4. We do not see any principal difficulties arising here, but it would take too much space to write the complete proof, so we restrict ourselves to this little observation only.

## 5 Volume of lattice simplices

For simplices we have a better method (than applying (8)) for bounding volume which appears in \[2, Theorem 3.4\] (see also \[3, Lemma 2.3\]). Let us reproduce this simple argument here.

###### Lemma 5

Let $`S=\text{conv}\{𝐯_0,\mathrm{},𝐯_d\}`$ be any simplex and let $`𝐰I_l(S)`$ have barycentric coordinates $`(\alpha _0,\mathrm{},\alpha _d)`$. Then

$$\text{vol}(S)\frac{l^d}{d!\times \alpha _1\times \alpha _2\times \mathrm{}\times \alpha _d}|I_l(S)|.$$
(28)
Proof. Consider the region

$$X=\{𝐰+_{i=1}^d\beta _i(𝐯_i𝐯_0)|\beta _i|\alpha _i1id\}.$$

It is a centrally symmetric parallelepiped around the vertex $`𝐰l^d`$ with volume $`\text{vol}(X)=d!\text{vol}(S)_{i=1}^d(2\alpha _i)`$. We have to show that the volume of $`X`$ cannot exceed $`(2l)^d|I_l(S)|`$. If this is not true, then $`X`$ contains (besides $`𝐰`$) at least $`|I_l(S)|`$ pairs of vertices $`𝐰\pm 𝐮l^d`$ by Corput’s theorem  and, clearly, at least one vertex of each such pair lies within $`I_l(S)`$, which is a contradiction.

Now we can deduce the following result.

###### Theorem 6

For any $`k1`$, the inequality (10) holds.

Proof. Let $`S^n`$ be lattice simplex with $`|I_l(S)|=k`$. By Theorem 2 there is $`𝐰I_l(S)`$ with $`m_S(𝐰)\gamma =(8l+7)^{2^{d+1}}/8`$. Assume $`\alpha _0\mathrm{}\alpha _d`$. Then it is easy to see that $`_{i=1}^d\alpha _i\gamma ^{d1}(1\gamma d)`$. The claim now follows from (28).

## 6 $`s(d,k,l)`$ for small $`d`$ and $`l`$

As we have already mentioned, $`s(2,k,1)`$ was computed by Scott . The simplex $`S_{2,k,l}`$ shows that $`s(2,k,l)l(l+1)^2(k+1)/2`$. Upper bounds on $`s(2,k,l)`$ can be obtained by applying (28) to the lower bounds on $`\beta (2,l)`$ from Table 1. But even if we knew $`\beta (2,l)`$ exactly, the best upper bound on $`s(2,k,l)`$ that this method would give is $`l^5k/2+O(l^4)k`$, so there would still be an uncertainty about $`s(2,k,l)`$.

Also, an interesting problem is the determination of $`s(3,k,1)`$. The simplex $`S_{3,k,1}`$ shows that $`s(3,k,1)6(k+1)`$. Theorem 6 gives, already for such small $`d`$, very bad bounds. However, there is a very simple argument, following the lines of Section 3 and proving that

$$s(3,k,1)\frac{29791}{2112}k<14.106k.$$
(29)
Given a lattice simplex $`S^3`$, we can deduce as before that the barycentric coordinates $`(\alpha _0,\mathrm{},\alpha _3)`$ of a lattice vertex maximizing $`m_S`$ satisfy $`2\alpha _31\alpha _0`$, $`3\alpha _21\alpha _0`$ and either $`4\alpha _21\alpha _0`$ or $`4\alpha _32\alpha _0`$. These inequalities do not guarantee yet that $`\alpha _0>0`$, but they guarantee that $`\alpha _1\frac{2}{31}`$ and, as it is routine to see, that $`\alpha _1\alpha _2\alpha _3\frac{2}{31}\times \frac{11}{31}\times \frac{16}{31}`$, which implies (29) by (28).

Of course, we could write more equations on the $`\alpha `$’s, but this method would not lead to the best possible bound. For example, the simplex $`S_{3,1,1}`$ shows that we cannot guarantee a vertex in $`I_1(S)`$ with $`\alpha _1\alpha _2\alpha _3>\frac{1}{2}\times \frac{1}{3}\times \frac{1}{12}`$, so the best bound we would hope to obtain this way is $`s(3,k,1)12k`$ only.

## 7 Lawrence’s Finiteness Theorem

A result of Lawrence \[4, Lemma 5\] implies that there is $`\gamma =\gamma (d)>0`$ such that, for any lattice simplex $`S=\text{conv}\{𝐯_0,\mathrm{},𝐯_d\}`$, the set $`I_1(S)`$ (if non-empty) contains a vertex $`𝐰`$ with $`\alpha _0\gamma `$, where $`(\alpha _0,\mathrm{},\alpha _d)`$ are the barycentric coordinates of $`𝐰`$. This directly follows from our Theorem 2 but unfortunately I could not find a simple argument giving the converse implication. And, in fact, the corresponding extremal functions are different: for example, $`\beta (2,1)1/8`$ while it is claimed in \[4, p. 439\] that we can take $`\gamma (2)=1/6`$.

However, one can deduce from  that $`\beta (d,1)>0`$ using the following simple modification of Lawrence’s proof. For $`i`$ let $`U_i=\{𝐰\mathrm{\Delta }_d\text{ca}(𝐰,\mathrm{\Delta }_d)<i\}`$, where $`\mathrm{\Delta }_d=\text{conv}\{\mathrm{𝟎},𝐞_1,\mathrm{},𝐞_d\}`$. Clearly, $`_{i=d+1}^{\mathrm{}}U_i=\text{int}(\mathrm{\Delta }_d)`$. By \[4, Theorem 3\] there exists $`i`$ such that for any $`𝐰\text{int}(\mathrm{\Delta }_d)`$ there are $`j`$ and $`𝐮^d`$ with with $`j𝐰+𝐮U_i`$. We claim that $`\beta (d,1)\frac{1}{i+1}`$. Indeed, let $`S^d`$ be any lattice simplex and let $`𝐯I_1(S)`$. Choose any affine function $`f:^d^d`$ with $`f(S)=\mathrm{\Delta }_d`$ and let $`𝐰=f(𝐯)\text{int}(\mathrm{\Delta }_d)`$. Given $`𝐰`$, choose the corresponding $`j`$ and $`𝐮^d`$. It is easy to check that $`𝐯^{}=f^1(j𝐰+𝐮)`$ belongs to $`I_1(S)`$ and satisfies $`m_S(𝐯^{})\frac{1}{i+1}`$.

But as it has already been remarked, Lawrence’s argument does not give any explicit bound.

## Acknowledgments

I am grateful to Günter Ziegler for drawing my attention to the problem and for very useful discussions and suggestions at all stages of this research. Also, I want to thank Marc Pfetsch and Julian Pfleifle for help with programming and the whole ”Discrete Geometry” group of TU-Berlin for their hospitality.

## C source code

```

/*

program for proving lower bounds on beta(d,l), the function defined in
e-print O.Pikhurko "Lattice Points inside Lattice Polytopes" at arXiv.org

Copyright (C) 2000 Oleg Pikhurko

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details:
http://www.gnu.org/copyleft/gpl.html

Purpose

For a lattice simplex $S\subset R^d$ and $x\in S$, let $m(x)$ be the
value of the smallest of the $d+1$ barycentric coordinates of
$x$. There is $\beta(d,l)>0$ such that if $S$ contains at least one
vertex of $lZ^d$ in its interior, then it contains such a vertex $x$
with $m(x)\ge \beta(d,l)$. The program computes lower bounds on
$\beta(d,l)$ for small $d$ and $l$.

Algorithm

Main idea: let $x\in S$ be a pair (almost) attaining $\beta(d,l)$. Let
$(a_0\le\dots\le a_d)$ be the barycentric coordinates of $x$. For any
numbers $p_0,\dots,p_d\in N$, there is $i\in[0,d]$ such that
$(lq+1)a_i-lp_i\le a_0$, where $q=\sum_{i=0}^d p_i$. (For otherwise we
have a contradiction, the vertex $x’=(q+1)x-\sum_{i=1}^d p_i v_i\in
lZ^d$ belongs to the interior of $S$ and has larger $m$, where $v_i$’s
are the vertices of $S$.)

That is, the $a_i$’s (for extremal $x$) must satisfy certain
linear inequalities. If these imply a lower bound on $a_0$ it is also
a lower bound on $\beta(d,l)$.

In the algorithm, we define our initial LP to consists of $0\le
a_0\le\dots\le a_d$ and $\sum_{i=0}^d a_i=1$. Then we repeat the
following. Given LP, let $(a_i)$ be a solution minimizing $a_0$. Try
to find $p_i$ with $(lq+1)a_i-lp_i>a_0$ which show that $(a_i)$ cannot
be the coordinates of extremal $x$, add the corresponding constraints
to our LP, and repeat.

Of all $p_i$’s, we choose the smallest in the lex order with minimum
$q$, which probably speeds convergence. Any better ideas?

As each time we have or-connected inequalities, we introduce binary
variables; of course, if $p_i=p_{i-1}$ then there is not need to
include the $i$th inequality. Also the $0$th equality is never
included because $(a_i)$ minimizes $a_0$ given LP, hence if
$(lq+1)a_0-lp_0>a_0$ for this LP, it will be true for any larger LP.

Usage

Prepare file "d-l-in.lp(s)" containing the lp program at which you
wish to start. (If the file does not exist, the program creates the
default initial lp.) Tip: rename file "d-l-last.lp(s)" as
"d-l-in.lp(s)" to start at the place, where the program finished last
time.

Run the program. Enter d l when asked. As each subsequent prompt,
enter 0 for non-iteractive mode, -1 to finish calculations, -2 to save
lp to the file "d-l-out.lp(s)" or a positive number which gives the
number of iterations to do till next prompt. At each termination of
the program (except by ^C), the last LP problem is saved to file
"d-l-last.lp(s)".

The program can be linked with either cplex (www.cplex.com) or lp_solve
(ftp://ftp.ics.ele.tue.nl/pub/lp_solve/). With this problem lp_solve is
unfortunately not very stable. Required libraries:
  cplex, socket, m, nsl (if compiled with -DCPLEX);
  liblpk.a, m, l (if compiled with the -DLPS).

Bugs

If we compile with -DLPS and read the initial lp from file, then
change_lp causes segmentation fault. However, everything works fine if
we compile with -DCPLEX or the initial problem is created by make_lp().

*/

#include <stdio.h>
#include <math.h>
#include <string.h>

/* library specific declarations; some of our function also depend
   on the library; these are placed at the end of the file */
#ifdef CPLEX

 #include <cplex.h>
 #define REAL double
 REAL tolerance=(REAL)1e-9; /* min tolerance allowed by CPLEX */
 REAL delta=(REAL)1e-9; /* initial upper bound on possible error on a[0];
                    for CPLEX it will be recomputed with each iteration */

 /* extension of all lp files */
 #define EXT ".lp"
 CPXENVptr env = NULL;
 CPXLPptr lp = NULL;
 int status = 0;

#elif defined LPS

 #include <lpkit.h>

 /* set different extension to avoid interference */
 #define EXT ".lps"
 lprec *lp;
 REAL delta=(REAL)1e-5; /* this is our guess for lp_solve */

#else

 #error Please compile with either -DCPLEX or -DLPS option

#endif

#define DMAX 9  /* size allocated for a[];
                   d must not exceed DMAX */
#define QMAX 1000 /* if next_p we get as far as q>=QMAX,
                     we assume that no approximation exists */
#define IMAX 600 /* the maximal number of iterations */

/* most important functions */
void get_a(void), read_a(void); /* get/(read from stdin) next solution a[] */
void next_p(void); /* compute approximation p[] */
void change_lp(void); /* add the corresponding constraints to lp */

/* auxiliary functions */
void print_a(REAL*), print_p(int*), print_all(void); /* various printing */
void die(char*); /* terminate the program */
void close_all(void); /* close all files, etc */
void init_lp(void); /* initialization */
void add_dl(char*);
void find_jumps(void);
void p_string(char*,int);
void save_lp(char*);
void check_solution(void);
void change_lp_plain(void);
void change_lp_tricky(void);

int d=2,l=1;
REAL best=(REAL) 0.0; /* the best lower bound we can guarantee so far */
REAL e=(REAL) 1e-6; /* small constant to be used in crucial roundings,
                    we search for p[] with (l*q+1)a[i]-l*p[i] > a[0]+e  */
REAL a[DMAX+1], b[DMAX+1]; /* b[] = (l*q+1)a[]-l*p[] */
int p[DMAX+1]; /* a, b & p have d+1 elements each */
int q;  /* q=p[0]+...+p[d] */
int jp, jp_ind[DMAX+1]; /* the number and array of jumps in p[] */
int it; /* the number of iterations */
int  d_stop; /* gives number of iteration before next prompt;
                if <=0 used for various actions */

int main()
{

  int next_stop=0;

  printf("Please enter d l ");
    scanf("%d %d",&d,&l);

  if (d>DMAX || d<=0 || l<=0)
    die("Wrong input values for d and l.");

  init_lp();

  for(it=0;it<IMAX;it++) {

    get_a(); /* read_a(); */

    check_solution(); /* can be commented out */

    if(best<a[0]-delta) /* we can improve our lower bound! */
      best=a[0]-delta;

    next_p();

    print_all();

    change_lp();

    if(it>=next_stop) {

      printf("\n\nEnter 0 (non-stop), -1 (end) -2 (save lp)\n");
      printf("or number (>0) runs till stop ");
      scanf("%d",&d_stop);

      if(d_stop==-1)
        die("Exiting...");

      if(d_stop==-2) {
        char lp_name[]="x-x-out" EXT;
        save_lp(lp_name);
      }

      if(d_stop==0)
        next_stop=IMAX+1;

      if(d_stop>0)
        next_stop=it+d_stop;

    }

  }

  die("The maximum number of iterations is exceeded.");

  return 0;

}

/* reads a[] from stdout; for testing purposes */
void read_a(void)
{

  int i;
  REAL suma;

  printf("Enter a\n");

  for(suma=0,i=0;i<d;i++) {
    scanf("%lf", &a[i]);
    suma+=a[i];
  }

  a[d]=1-suma;

  if (a[0]<0 || a[d]<0)
    die("a[] is not feasible.");

}

/* for given a[] compute a smallest p[] such that min((lq+1)a[]-l*p[])>a[0].
   return 0 if no reasonably small P exists */
void next_p(void)
{

  int found=0,i,sump,diff;

  /*  printf("next_p: obtained the following vector:\n");
  print_a(a); */

  q=0; /* we try to find p[] with smallest q, so we interate on q */

  do {

    q++;

    sump=0;
    for(i=0;i<=d;i++) {

      /* the case i==0 receives a special treatment */
      p[i]=ceil(q*(a[i]-(i?0:delta)) + (a[i]-a[0])/l - 1 - e);

      if (p[i]<0) /* this may happen if a[i] is small */
        p[i]=0;
      sump+=p[i];
    }

    if (sump>=q) {

      found=1;

      if (sump>q) {  /* here different p[] are possible;
                        we take the lex-smallest p[] */
        diff=sump-q;
        for(i=0; i<=d && diff>0; i++)
          if(p[i]<=diff) {
            diff-=p[i];
            p[i]=0;
          } else {
            p[i]-=diff;
            diff=0;
          }
      }

    }

  } while (! found && q< QMAX);

  if (found){

    /* let us run small check; just in case */
    for(sump=0,i=0;i<=d;i++)
      sump+=p[i];

    if (sump != q)
      die("next_p: sump does not equal q");

    for(i=0;i<=d;i++) {
      b[i]=(l*q+1)*a[i]-l*p[i];
      if( p[i]<0 || ( b[i]<= a[0] && a[i]>0 ) )
        die("next_p: the sequence p[] is bad");
    }

  } else
    die("Suitable p cannot be found. Success!?");

}

/* creates a string "r_p[0]_..._p[d]_i" */
void p_string(char s[], int i)
{

  int j;
  char temp_name[255];

    strcpy(s,"r");
    for(j=0;j<=d;j++) {
      sprintf(temp_name,"_%d",p[j]);
      strcat(s,temp_name);
    }
    sprintf(temp_name,"_%d",i);
    strcat(s,temp_name);

}

/* modifies our lp */
void change_lp(void)
{

  find_jumps(); /* compute the jump in p[] */

  if(!jp) /* a check just in case */
    die("No jumps in p[], which cannot happen at all.");

  change_lp_tricky();

  /* change_lp_plain(); */

  /* change_lp_plain() produces lp which are longer to solve, but
     whose binary variables are easier to understand (by humans);
     we keep the definition of change_lp_plain() for debugging */

}

/* compute the number jp and array jp_ind[] of jumps in p[]  */
void find_jumps(void)
{

  int i,p_prev;

  jp=0;
  p_prev=p[0];

  for(i=1;i<=d;i++)
    if (p[i]>p_prev) {
      jp_ind[jp++]=i;
      p_prev=p[i];
    }

}

/* prints d+1 array of reals */
void print_a(REAL s[])
{

  int i;

  for(i=0;i<=d;i++)
    printf("%19.17f ",s[i]);
  printf("\n");

}

/* prints d+1 array of integers */
void print_p(int s[])
{

  int i;

  for(i=0;i<=d;i++)
    printf("%d ",s[i]);
  printf("\n");

}

/* prints some global information */
void print_all()
{

  printf("\nGlobal variables:\n");
  printf("d=%d, l=%d, iteration=%d, q=%d, delta=%.2e\n",
         d,l,it,q,delta);

#ifdef CPLEX
  if(status) { /* get the error message corresponding to status */
      char  errmsg[1024];
      CPXgeterrorstring (env, status, errmsg);
      printf ("%s\n", errmsg);
  }
#endif

  printf("Best guaranteed lower bound so far = %19.17f\n", best);
  printf("Vectors a[] p[] & b[]:\n");
  print_a(a);
  print_p(p);
  print_a(b);

}

/* stop execution */
void die(char s[])
{

  printf("\n%s\n",s);
  print_all();

  close_all();

  exit(1);

}

/* replaces the 1st and 3d latters in a filename by the values of d and l */

void add_dl(char s[])
{

  s[0]=’0’+d;
  s[2]=’0’+l;

}

/* check if new a[] is compatible with p[]; lp_solve sometimes returns
   an unfeasible solution and check_solution fails - any ideas? */
void check_solution(void)
{

  int i,found=0;

  for(i=0;i<=d;i++)
    if((b[i]=(l*q+1)*a[i]-l*p[i]) <= a[0]+e)
      found=1;

  if (!found)
    die("check_solution failed.");

}

/* library specific definitions of functions: change_lp_plain(),
                change_lp_tricky(), close_all(), save_lp() */

#ifdef CPLEX

/* Initialization */
void init_lp(void)
{

  int i;

  char in_lp[]="x-x-in" EXT;
  FILE *lp_file;

  add_dl(in_lp);

  /* Initialize the CPLEX environment */
  env=CPXopenCPLEXdevelop(&status);
  if (env == NULL)
    die("CPXopenCPLEXdevelop failed.");

  /* sets extra space for columns */
  status=CPXsetintparam(env,CPX_PARAM_COLGROWTH,DMAX*IMAX);
  if(status)
    die("Cannot set CPX_PARAM_COLGROWTH.");

  /* sets extra space for rows */
  status=CPXsetintparam(env,CPX_PARAM_ROWGROWTH,(DMAX+1)*IMAX);
  if(status)
    die("Cannot set CPX_PARAM_ROWGROWTH.");

  /* sets extra space for non-zero elements;
     do we need this, given two settings above? */
  status=CPXsetintparam(env,CPX_PARAM_NZGROWTH,(DMAX+1)*(DMAX+1)*IMAX);
  if(status)
    die("Cannot set CPX_PARAM_NZGROWTH.");

  /* sets integrality tolerance to smallest allowable */
  status=CPXsetdblparam(env,CPX_PARAM_EPINT,tolerance);
  if(status)
    die("Cannot set CPX_PARAM_EPINT.");

  /* sets optimallity tolerance to smallest allowable */
  status=CPXsetdblparam(env,CPX_PARAM_EPOPT,tolerance);
  if(status)
    die("Cannot set CPX_PARAM_EPOPT.");

  /* sets feasibility tolerance to smallest allowable */
  status=CPXsetdblparam(env,CPX_PARAM_EPRHS,tolerance);
  if(status)
    die("Cannot set CPX_PARAM_EPRHS.");

  /* create lp environment */
  lp=CPXcreateprob(env,&status,"plattice");
  if (lp == NULL)
    die("CPXcreateprob failed.");

  if((lp_file=fopen(in_lp, "r"))==NULL) {

    REAL obj[DMAX+1]; /* objective function */
    char colname[DMAX+1][4]; /* names of variables, "a0", "a1", etc */
    char* colname_p[DMAX+1];
    REAL rhs[DMAX+1],rmatval[3*DMAX+1];
    int rmatbeg[DMAX+2],rmatind[3*DMAX+1];
    char sense[DMAX+1];

    printf("File %s does not exists.\n",in_lp);
    printf("Creating the basic initial lp.\n");

    for(i=0;i<=d;i++) {
      obj[i]=(REAL)0.0;
      colname_p[i]=strcpy(colname[i],"a0");
      colname[i][1]+=i;
    }

    obj[0]=(REAL)1.0;

    /* add real, >=0 (by default) variables a0,a1,... */
    status=CPXnewcols(env,lp,d+1,obj,NULL,NULL,NULL,colname_p);
    if (status)
      die("CPXnewcols failed.");

    for(i=0;i<d;i++) { /* inequality -a[i] + a[i+1] >= 0 */
      rhs[i]=(REAL)0.0;
      sense[i]=’G’;
      rmatbeg[i]=2*i;
      rmatind[2*i]=i;
      rmatval[2*i]=(REAL)-1.0;
      rmatind[2*i+1]=i+1;
      rmatval[2*i+1]=(REAL)1.0;
    }

    /* equation a0+a1+...=1 */
    rhs[d]=(REAL)1.0;
    sense[d]=’E’;
    rmatbeg[d]=2*d;
    for(i=0;i<=d;i++) {
      rmatind[2*d+i]=i;
      rmatval[2*d+i]=(REAL)1.0;
    }

    /* add these basic constraints */
    status=CPXaddrows(env,lp,0,d+1,3*d+1,rhs,sense,
                      rmatbeg,rmatind,rmatval,NULL,NULL);
    if (status)
      die("CPXaddrows failed.");

    } else {

    fclose(lp_file);

    status=CPXreadcopyprob(env,lp,in_lp, NULL);
    if (status)
      die("CPXreadcopyprob failed.");

  }

}

/* compute a[] as the next solution; */
void get_a(void)
{

  /* solve the lp */

  status=CPXmipopt(env, lp);
  if (!status) {

    status = CPXgetmipx(env, lp, a, 0, d);
    if (status)
      die("CPXgetmipx failed\n");

  } else { /* CPXmipopt failed; maybe there are no bin variables yet */

    /*    printf("%d\n",status); */

    status=CPXprimopt(env,lp);
    if(status)
      die("Both CPXmipopt and CPXprimopt failed.");
    status=CPXgetx(env,lp,a,0,d);
    if(status)
      die("CPXgetx failed.");

  }

}

/* add constraints corresponding to p[] to the lp;
   a straightforward version */
void change_lp_plain(void)
{

  int j,i;
  int cur_numcols; /* the number of columns (variables) in lp */
  /* variables for CPXaddcols */

  REAL ub[DMAX];

  /* variables for CPXaddrows;
     we may need as much as d+1 equations */
  REAL rhs[DMAX+1];
  int rmatbeg[DMAX+2];
  int rmatind[4*DMAX];
  char sense[DMAX]; /* types of constraints */
  char row_name[DMAX+1][255];
  char* row_pointer[DMAX]; /* symbolic names of rows */
  REAL rmatval[4*DMAX]; /* our equations have <=4*d non-zero coefficients */
  char ctype[DMAX]; /* types of variables */

  /* first we add jp binary variables z[0],...,z[jp-1] */
  /* for j=1 we add an unecessary binary variable,
     but it is fixed to 1 anyway and brings no harm;
     we must also specify ub[j]=1 because otherwise
     CPXmipopt returns 3002. */
  for(j=0;j<jp;j++) {
    ctype[j]=’B’;
    ub[j]=(REAL)1.0;
  }
  status=CPXnewcols(env,lp,jp,NULL,NULL,ub,ctype,NULL);
  if(status)
    die("CPXnewcols failed.");

  /* get the number of columns (variables) in lp */
  cur_numcols = CPXgetnumcols (env, lp);
  if(cur_numcols<=0)
    die("CPXgetnumcols failed.");

  /* now out added variables have numbers
     from cur_numcols-jp to cur_numcols-1 */

  /* now we add jp constraints of the form
     -a[0]+(l*q+1)*a[i]+(l*q+1) z[j]<= l*p[i]+(l*q+1), where i=jp_ind[j] */

  if(delta<tolerance*(3+2*l*q)) /* recompute the max abs mistake */
    delta=tolerance*(3+2*l*q);

  for(j=0;j<jp;j++) {

    i=jp_ind[j];

    rhs[j]=(REAL) l*p[i]+(l*q+1);
    rmatbeg[j]=3*j; /* three non-zero coeff. per row */
    rmatval[3*j]=(REAL) -1;
    rmatind[3*j]=0; /* x[0] */
    rmatval[3*j+1]=(REAL) l*q+1;
    rmatind[3*j+1]=i; /* x[i] */
    rmatval[3*j+2]=(REAL) l*q+1;
    rmatind[3*j+2]=cur_numcols-jp+j; /* z[j] */
    sense[j]=’L’;
    p_string(row_name[j],i);
    row_pointer[j]=row_name[j];

  }

  /* plus the constraint z[0]+...+z[jp-1]>=1; */

  if(delta<tolerance*jp) /* recompute the max abs mistake */
    delta=tolerance*jp;

  rhs[jp]=(REAL)1.0;
  sense[jp]=’G’;
  rmatbeg[jp]=3*jp;
  rmatbeg[jp+1]=4*jp;
  p_string(row_name[jp],0);
  row_pointer[j]=row_name[j];

  for(j=0;j<jp;j++) {
    rmatval[3*jp+j]=(REAL) 1;
    rmatind[3*jp+j]=cur_numcols-jp+j;
  }

  status=CPXaddrows(env,lp,0,jp+1,4*jp,rhs,sense,
                    rmatbeg,rmatind,rmatval,NULL,row_pointer);
  if(status)
    die("CPXaddrows failed.");

}

/* add constraints corresponding to p[] to the lp; uses the minimum
   nnumber of binary variables; a bit tricky */
void change_lp_tricky(void)
{

  /* variables for CPXaddrows;
     we may need as much as d equations */
  REAL rhs[DMAX];
  int rmatbeg[DMAX+1];
  int rmatind[DMAX*(DMAX+2)]; /* non-zero elements in added equations */
  char sense[DMAX]; /* types of constraints */
  char row_name[DMAX][255];
  char* row_pointer[DMAX]; /* symbolic names of rows */
  REAL rmatval[DMAX*(DMAX+2)];

  int j, jp_log=0, j2=1;
  int cur_numcols; /* the number of columns (variables) in lp */

  /* we add constraints saying that there is j<jp such that
     (l*q+1)a[i]- a[0]<= l*p[i], where i=jp_ind[j];
     to do so we introduce jp_log binary variables, where
     jp_log is the smallest number with $2^jp_log>= jp.

     For each of 2^jp_log possible combinations of binaries all
     inequalities are vacuous except one, which is precisely one of
     our desired inequalities. We do not explain how such a system
     is constructed but encourage the reader to check the algorithm
     and to run the program and look at the resulting lp */

  while(j2 < jp) {
    j2*=2;
    jp_log++;
  }

  if(delta<tolerance*(l*q+1)*(jp_log+1))  /* recompute the max mistake */
    delta=tolerance*(l*q+1)*(jp_log+1);

  if (jp_log>0) {

    /* variables for CPXaddcols */
    REAL ub[DMAX];
    char ctype[DMAX]; /* types of variables */

    /* first we add jp_log binary variables
       we must also specify ub[j]=1 because otherwise
       CPXmipopt returns 3002. */

    for(j=0;j<jp_log;j++) {
      ctype[j]=’B’;
      ub[j]=(REAL)1.0;
    }

    status=CPXnewcols(env,lp,jp_log,NULL,NULL,ub,ctype,NULL);
    if(status)
      die("CPXnewcols failed.");

  }

  /* now we add constraints, the case jp_log==0 inclusive */

  /* get the number of columns (variables) in lp */
  cur_numcols = CPXgetnumcols (env, lp);
  if(cur_numcols<=0)
    die("CPXgetnumcols failed.");
  /* our added variables have numbers
     from cur_numcols-jp_log to cur_numcols-1 */

  for(j=0;j<jp;j++) {

    int jj,jj2,i;

    i=jp_ind[j];

    sense[j]=’L’;

    rmatbeg[j]=j*(2+jp_log); /* 2+jp_log non-zero entries per row */

    rmatind[j*(2+jp_log)]=0; /* coefficient at a0 */
    rmatval[j*(2+jp_log)]=(REAL)-1.0;

    rmatind[j*(2+jp_log)+1]=i; /* coefficient at a[i] */
    rmatval[j*(2+jp_log)+1]=(REAL) (l*q+1);

    p_string(row_name[j],i);
    row_pointer[j]=row_name[j];

    rhs[j]=(REAL) l*p[i];

    for(jj=0,jj2=1; jj<jp_log; jj++,jj2*=2) {

      /* coefficient at jj’th binary variable */
      rmatind[j*(2+jp_log)+2+jj]=cur_numcols-jp_log+jj;

      if( ( (unsigned int) 1 << jj) & j ) {

        rhs[j]+=(REAL)l*q+1;
        rmatval[j*(2+jp_log)+2+jj]=(REAL)+l*q+1;
      } else {
        if (j+jj2>=jp)
          rmatval[j*(2+jp_log)+2+jj]=(REAL)0.0;
        else
          rmatval[j*(2+jp_log)+2+jj]=(REAL)-l*q-1;
      }
    }

  }

  status=CPXaddrows(env,lp,0,jp,jp*(2+jp_log),rhs,sense,
                    rmatbeg,rmatind,rmatval,NULL,row_pointer);
  if(status)
    die("CPXaddrows failed.");

}

/* close log-files, cplex and write the current program to file */

void close_all(void)
{

  char out_lp[]="x-x-last" EXT;

  save_lp(out_lp);
  CPXwriteprob(env,lp,out_lp, NULL);
  CPXfreeprob(env, &lp);
  CPXcloseCPLEX(&env);

}

/* saves the lp to file; we cannot here call die()
   because save_lp may be called from it */

void save_lp(char s[])
{

  add_dl(s);
  printf("Writing current LP to file %s ... ",s);
  status=CPXwriteprob(env,lp,s, NULL);

  if(status)
    printf("failed, status=%d.\n", status);
  else
    printf("ok.\n");

}

#elif defined LPS

/* read/create the initial lp */
void init_lp(void)
{

  nstring col_name;
  char in_lp[]="x-x-in" EXT;
  FILE *lp_file;

  add_dl(in_lp);

  if((lp_file=fopen(in_lp, "r"))==NULL) {

    int i,j;
    REAL row[DMAX+2];

    printf("File %s does not exists.\n",in_lp);
    printf("Creating the basic initial lp.\n");

    lp=make_lp(0,d+1);
    if(lp==NULL)
      die("make_lp failed.");

    for(i=0;i<=d;i++) { /* give variables names a0,a2,... */
      strcpy(col_name,"a0");
      col_name[1]+=i;
      set_col_name(lp,i+1,col_name);
    }

    /* constraint -a[i-1]+a[i]>=0 */
    /* add_constraint starts indexing with 1 and sometimes
       the 0th entry is also used for other purposes; so we
       have to modify all indexing accordingly */

    for(i=1;i<=d;i++) {
      for(j=0;j<=d+1;j++)
        row[j]=(REAL) 0.0;
      row[i]=(REAL) -1.0;
      row[i+1]=(REAL) 1.0;
      add_constraint(lp,row,GE,(REAL) 0.0);
    }

    /* a[0]+...+a[d]=1 */
    for(j=1;j<=d+1;j++)
      row[j]=(REAL) 1.0;
    add_constraint(lp,row,EQ, (REAL) 1.0);

    /* set objective function */

    row[1]=(REAL) 1.0;
    for(j=2;j<=d+1;j++)
      row[j]=(REAL) 0.0;

    set_obj_fn(lp,row);

    print_lp(lp);

  } else {

    /* NB: keep lp_file open */
    lp=read_lp_file(lp_file,FALSE,"plattice");
    fclose(lp_file);

    if(lp==NULL)
      die("read_lp_file failed.");

  }

}

/* compute a[] as the first d+1 variables of the solution; */
void get_a(void)
{

  int i;

  if(solve(lp))
    die("solve failed.");

  for(i=0;i<=d;i++)
    a[i]=lp->best_solution[lp->rows+1+i];

}

/* add constraints corresponding to p[] to the lp;
   a straightforward version */
void change_lp_plain(void)
{

  nstring row_name; /* ! lp_solve allows at most 24 characters */

  int j,i,jj;

  REAL col[(DMAX+1)*(IMAX+1)];
  /* (DMAX+1)*(IMAX+1) is max possible number of rows/cols in lp */

  for(j=0;j<=lp->rows;j++)
    col[j]=(REAL) 0.0;

  /* first we add jp binary variables z[0],...,z[jp-1] */
  /* for j=1 we add an unecessary binary variable,
     but it is fixed to 1 anyway and brings no harm; */

  for(j=0;j<jp;j++) {
    add_column(lp,col);
    set_int(lp,lp->columns,TRUE);  /* the variable is integer */
    set_upbo(lp,lp->columns,(REAL) 1.0);
    set_lowbo(lp,lp->columns,(REAL) 0.0);
  }

  /* now we add jp constraints of the form
     -a[0]+(l*q+1)a[i]+(l*q+1) z[j]<= l*p[i]+(l*q+1), where i=jp_ind[j] */
  for(j=0;j<jp;j++) {
    i=jp_ind[j];
    for(jj=0;jj<=lp->columns;jj++)
      col[jj]=(REAL) 0.0;
    col[1]=(REAL) -1.0; /* at a[0] */
    col[i+1]=(REAL) l*q+1; /* at a[i] */
    col[lp->columns-jp+j+1]= (REAL) l*q+1; /* at z[j] */
    add_constraint(lp,col,LE,(REAL) l*p[i]+(l*q+1));
    p_string(row_name,i);
    set_row_name(lp,lp->rows,row_name);
  }

  /* plus the constraint z[0]+...+z[jp-1]>=1 */
  for(jj=0;jj<=lp->columns;jj++)
    col[jj]=(REAL) 0.0;
  for(j=1;j<=jp;j++)
    col[lp->columns-jp+j]=(REAL) 1.0;
  add_constraint(lp,col,GE,(REAL) 1.0);
  p_string(row_name,0);
  set_row_name(lp,lp->rows,row_name);

}

/* add constraints corresponding to p[] to the lp; uses the minimum
   nnumber of binary variables; a bit tricky */
void change_lp_tricky(void)
{

  int j, jp_log=0, j2=1;

  nstring row_name; /* ! lp_solve allows at most 24 characters */

  REAL col[(DMAX+1)*(IMAX+1)];
  /* (DMAX+1)*(IMAX+1) is max possible number of rows/cols in lp */

  for(j=0;j<=lp->rows;j++)
    col[j]=(REAL) 0.0;

  /* we add constraints saying that there is j<jp such that
     (l*q+1)a[i]- a[0]<= l*p[i], where i=jp_ind[j];
     to do so we introduce jp_log binary variables, where
     jp_log is the smallest number with $2^jp_log>= jp.

     For each of 2^jp_log possible combinations of binaries all
     inequalities are vacuous except one, which is precisely one of
     our desired inequalities. We do not explain how such a system
     is constructed but encourage the reader to check the algorithm
     and to run the program and look at the resulting lp */

  while(j2 < jp) {
    j2*=2;
    jp_log++;
  }

  /* first we add jp_log binary variables */
  for(j=0;j<jp_log;j++) {
    add_column(lp,col);
    set_int(lp,lp->columns,TRUE);  /* the variable is integer */
    set_upbo(lp,lp->columns,(REAL) 1.0);
    set_lowbo(lp,lp->columns,(REAL) 0.0);

  }

  /* our added variables have numbers
     from lp->columns-jp_log to lp->columns-1 */

  /* add constraints now */

  for(j=0;j<jp;j++) {

    int jj,jj2,i,jj_beg;
    REAL rhs;

    for(jj=0;jj<=lp->columns;jj++)
      col[jj]=(REAL) 0.0;

    i=jp_ind[j];

    col[1]=(REAL) -1.0; /* coefficient at a0 */
    col[i+1]=(REAL) l*q+1; /* coefficient at a[i] */

    rhs=(REAL) l*p[i]; /* rhs as it were without binary variables */

    for(jj=0,jj2=1; jj<jp_log; jj++,jj2*=2) {

      /* coefficient at jj’th binary variable */
      jj_beg=lp->columns-jp_log+jj+1;
      if( ( (unsigned int) 1 << jj) & j ) {
        rhs+=(REAL)l*q+1;
        col[jj_beg]=(REAL) l*q+1;
      } else {
        if (j+jj2>=jp)
          col[jj_beg]=(REAL) 0.0;
        else
          col[jj_beg]=(REAL) -l*q-1;
      }
    }

    add_constraint(lp, col, LE, rhs);
    p_string(row_name,j);
    set_row_name(lp,lp->rows,row_name);
  }

}

/* write out d-l-last.lp(s) */
void close_all(void)
{

  char out_lp[]="x-x-last" EXT;

  save_lp(out_lp);

}

/* saves the lp to file; we cannot here call die()
   because save_lp may be called from it */

void save_lp(char s[])
{

  FILE *lp_file;

  add_dl(s);
  printf("Writing current LP to file %s ... ",s);

  if((lp_file=fopen(s,"w"))==NULL)

    printf("failed.\n");

  else {

    write_LP(lp,lp_file); /* purify complains here about "uninitialized memory
                             read", though the code seems clean to me */

    if(fclose(lp_file))
      printf("failed.\n");
    else
      printf("ok\n");

  }

}

#endif

```


text:
# Fermion transmutation—a renormalization effect in gauge theory

## 1 Introduction

By ‘fermion transmutation’ we mean a process in which a fermion changes its generation index as a direct consequence of the rotation of the mass matrix, and not as a secondary effect such as $`e\mu `$ conversion via FCNC. Examples of transmutation are: $`e\mu ,e\tau ,\mu \tau `$, which can occur e.g. in a Compton-like process schematically represented in Figure 1.

That the fermion mass matrix rotates, by which we mean that it undergoes unitary transformations as the scale changes, can be seeen from the renormalization group equation, both in the standard model (SM) and in the dualized standard model (DSM).

An earlier talk summarizes our earlier work on DSM. In this talk I shall report on two effects, namely transmutational decays and photo-transmutations, in both SM and DSM, but with emphasis on DSM.

## 2 Mass matrix rotation

The SM renormalization group equation for the charged lepton mass matrix $`L`$ has a term which, given that the leptonic MNS mixing matrix $`U`$ is nontrivial, rotates it as the scale $`\mu `$ changes. The linearized equation:

$$\frac{dL}{d\mu }=\frac{3}{128\pi ^2}\frac{1}{246^2}(ULU^{})(ULU^{})^{}L+\mathrm{},$$
(1)
where $`ULU^{}=N`$ the neutrino (Dirac) mass matrix, already shows that even if $`L`$ is diagonal at a chosen scale it cannot remain so at all other scales. The magnitude of the off-diagonal elements will depend on poorly known or unknown quantities such as the mixing $`U`$ and the Dirac mass $`m_3`$ of the heaviest neutrino. If we take the present popular theoretical biases of $`U`$ bimaximal and $`m_3m_t`$, then (1) gives for each decade change in energy:

$`\mu |\tau `$ $`\mathrm{changes}\mathrm{by}`$ $`5.5\times 10^3\mathrm{GeV}`$
$`e|\tau `$ $`\mathrm{changes}\mathrm{by}`$ $`1.8\times 10^7\mathrm{GeV}`$
$`e|\mu `$ $`\mathrm{changes}\mathrm{by}`$ $`1.1\times 10^8\mathrm{GeV}`$
In the DSM, the fermion mass matrix is of the following factorized form:

$$m=m_T\left(\begin{array}{c}x\\ y\\ z\end{array}\right)(x,y,z),$$
(2)
where $`m_T`$ is essentially the mass of the heaviest generation. Under renormalization $`m`$ remains factorized, but the vector $`(x,y,z)`$ changes as

$$\frac{d}{d\mu }\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\frac{5}{32\pi ^2}\rho ^2\left(\begin{array}{c}x_1\\ y_1\\ z_1\end{array}\right),$$
(3)
where $`\rho `$ is a (fitted) constant and

$$x_1=\frac{x(x^2y^2)}{x^2+y^2}+\frac{x(x^2z^2)}{x^2+z^2},\mathrm{cyclic}.$$
(4)
The off-diagonal elements have been calculated explicitly, using 3 free parameters determined by fitting experimental mass and mixing parameters (giving sensible predictions for the remaining paramenters). These are shown in Figure 2. Hence the results we report below are entirely parameter-free.

## 3 Definition of lepton states

To define lepton states one must diagonalize the mass matrix, but since the eigenvectors depend on scale, there is no canonical recipe for doing so.

We suggest two quite different schemes for exploration. In the fixed scale diagonalization scheme (FSD) the state vectors are fixed at a chosen scale. As a result, the mass matrix is diagonal only at that scale. Such a scheme is applicable to SM, where it is tacitly assumed in most calculations.

For the DSM we use a more sophisticated scheme which we may call step-by-step diagonalization (SSD). It is a working criterion in which the lepton state vectors are always orthogonal and the mixing matrix is always unitary. There are 3 main steps:
$``$ run the $`3\times 3`$ mass matrix to a scale equal to mass of heaviest generation;
$``$ corresponding eigenvector is its state vector;
$``$ repeat with the remaining $`2\times 2`$ submatrix.
In contrast to FSD, the charged lepton mass matrix here is diagonal at scales of $`m_\tau ,m_\mu ,m_e`$, as seen in Figure 2.

## 4 Transmutational decays

With the lepton states defined as in §3 and the rotations obtained in §2 we can now study transmutational processes, the most obvious of which are decays.

By expanding the fermion propagator we get for $`\mathrm{}_\alpha \mathrm{}_\beta ,\alpha \beta `$:

$$\mathrm{transmutation}/\mathrm{diagonal}\alpha |\beta /E,$$
(5)
where $`E`$ is a typical energy for the decay. Estimates are in Table 1.

With significant exceptions, SM (with FSD) estimates are all far below present experimental bounds and are hence not so interesting. An exception is the process $`\mu ^{}e^{}e^+e^{}`$, where one could get a branching ratio of $`10^3`$ (limit $`10^{12}`$), if one applied FSD naively.

The parameter-free calculations in DSM give branching ratios which are in general larger but still below present experimental limits. It is important to note that because of SSD the branching ratios of transmutational leptonic decays are automatically zero to first order (Figure 2). The $`\pi ^0`$ decay is of particular interest as being less than one order from the experimental bound.

## 5 Photo-transmutation

We have studied the following process, mainly for DSM:

$$\gamma +\mathrm{}_\alpha \gamma +\mathrm{}_\beta ,\alpha \beta .$$
(6)
There are two points to note. First, the kinematics is unfamiliar, with the mass eigenstates being varying linear combinations of the $`\tau ,\mu ,e`$ states. The other point is that standard formulae for summing $`\gamma `$ matrix traces cannot be applied directly. We calculated the individual spin and polarization amplitudes and then summed them by hand.

We calculated the cross sections for: $`\gamma e\gamma \mu ,\gamma e\gamma \tau ,\gamma \mu \gamma \tau `$, leaving out $`\tau `$-initiated reactions as experimentally unrealistic at present.

A sample of the DSM results is presented in Figure 3. We get in general:

$$\gamma \mu \gamma \tau >\gamma e\gamma \tau >\gamma e\gamma \mu .$$

However, at low energies $`\gamma e\gamma \mu `$ becomes quite sizeable, Figure 4, with the total cross section having a peak of $``$ 100 pb at c.m. energy $``$ 200 MeV.

Calculations in SM depend on further assumptions. If the Dirac masses of the neutrinos are hierarchical, then above $`\tau `$ the cross sections are approximately given by scaling those of DSM with the relevant rotation matrix elements. For example $`\gamma \mu \gamma \tau `$ at $`\sqrt{s}=17.8`$ GeV is $`23`$ orders smaller.

## 6 Possible experimental tests in the near future

Since only the DSM estimates and calculations are parameter-free, while the SM results presented here are subject to further assumptions and uncertainties, we shall point out tests for DSM only here.

The estimates for the transmutational decay modes: $`\pi ^0\mu ^{}e^+,\psi \tau ^{}\mu ^+,\mathrm{{\rm Y}}\tau ^{}\mu ^+`$ could be near experimental limits and sensitivites, for LEP, BEPC and B-factories.

For photo-transmutations, one may consider virtual $`\gamma `$ from $`e^+e^{}`$ colliders. Above $`\tau ,e\tau `$ is more important, while below it is $`e\mu `$. Again, LEP and/or BEPC may provide tests.

## 7 Conclusions

$``$ Fermion transmutation necessarily occur in both SM and DSM.
$``$ The SM results are in general smaller than DSM, with uncertainties.
$``$ The DSM calculations are entirely parameter-free. There are no violations of data in all the cases we were able to consider.
$``$ Experimental tests of DSM predictions seem feasible in the near future: (1) decays of $`\pi ^0,\psi ,\mathrm{{\rm Y}}`$, (2) photo-transmutation of $`\gamma e\gamma \tau `$ at high and $`\gamma e\gamma \mu `$ at low energy, (3) other processes e.g. $`e^+e^{}e^+\mu ^{}`$.
$``$ Transmutation leads to exciting new physics that has to be explored.


text:
# Quantum Interference Mechanism of Cooperative Optical Phenomena in Extended Media

## Abstract

In the quantum process of stimulated Raman scattering (SRS), a laser photon propagating in a resonance medium undergoes multifold conversions into a Stokes photon and back. The nontrivial “cooperative” behavior of the Stokes component of light transmitted through the medium is proven to be completely determined by the interference of scattering amplitudes in different sub-channels of the Stokes channel, which obviously combines all the sub-channels with an odd number of photon conversions. The theory of superfluorescence is then derived as the limiting case of the SRS theory.

It has been typical to treat cooperative phenomena of resonance quantum optics in an extended medium, such as stimulated Raman scattering (SRS)  and superfluorescence (SF) , by first describing the initial growth of the Stokes pulse in SRS (or, respectively, the radiation emitted in a SF process) by linear quantum theory and then the succeeding evolution by the classical nonlinear equations. It is assumed that once the Stokes wave is initialized by the quantum vacuum fluctuations and the Stokes pulse has grown to a classical value, then the classical equations become valid. Both the classical Maxwell-Bloch (MB) equations describing SF  and a set of the classical equations describing SRS  were proven  to be integrable. That allows one to apply the inverse scattering method  to study the time evolution of the SF radiation and Stokes pulses within a medium.

In this Letter, we propose and discuss in detail a microscopic mechanism of cooperative optical phenomena based on studies of a quantum version of the standard classical model of SRS.

Propagating in a partially excited resonance medium, an incident laser photon undergoes multifold conversions into a Stokes photon and back. The number of these conversions defines a sub-channel of the scattering process. A scattering amplitude in the elastic (Rayleigh) channel is obviously given by a sum of scattering amplitudes in sub-channels with an even number of conversions. In this channel, an incident laser photon leaves a medium preserving the frequency. Correspondingly, a scattering amplitude in the inelastic (Stokes) channel is given by a sum of the scattering amplitudes in sub-channels with an odd number of conversions. Therefore, a photon leaves a medium as a Stokes one.

Analyzing explicit expressions for the scattering amplitudes, we show that a nontrivial “cooperative” behavior of the Stokes component of the radiation transmitted through a medium is completely determined by the interference of the scattering amplitudes in different sub-channels of the inelastic channel of scattering. The interference of the scattering amplitudes in different sub-channel of the elastic channel results correspondingly in the “cooperative” behavior of the Rayleigh component of transmitted light.

A solution of the quantum model of SRS shows that there is neither “linear quantum stage” nor “succeeding nonlinear classical evolution” of the Raman process. The quantum interference mechanism describes the “cooperative” character of SRS in a unified manner for all incident photons. Moreover, the scattering amplitudes do not depend on any physical parameters of an incident laser pulse, such as its amplitude, duration, etc. They are determined only by the quantum state of a medium. Scattering in a medium, even in the elastic channel, an incident laser photon changes the quantum state of a medium. This results in changing the scattering amplitudes for all successive photons, explaining thus the “cooperative” behavior of the transmitted light.

When the number of excitations of a medium grows, a number of possible inelastic sub-channels also grows, and the probability of scattering in the Stokes channel increases. However, if a number of excited atoms exceeds a half of the total number of atoms in a medium, further excitation results in a decrease in the number of inelastic sub-channels. Therefore, the probability of scattering in the Stokes channel begins to fall, and vanishes asymptotically in the limit of a completely inverted medium. The intensity of the Stokes radiation coming out of a medium has a pulse (or soliton-like) shape.

Finally, we show that a quantum model of SF in an extended medium (or a quantum version of the classical MB model) does not require a special consideration. Its solution is derived as the limiting case of a solution of the SRS problem with an appropriate initial condition.

The “cooperative” behavior of the radiation emitted by an initially excited system of two-level atoms is determined by quantum interference of amplitudes in different sub-channels of emission. In the SF case, the sub-channels of the emission process are characterized by different numbers of possible reabsorption processes of a photon propagating in a partially excited medium.

Integrability of the quantum model of SRS allows one to derive an exact expression for the out-state of a scattering process . Unfortunately, a complex algebra of operators involved in this expression does not admit an analytical computation of physical observables of the transmitted radiation. The problem of a computation of physical observables is equivalent to the problem of correlation functions in the theory of integrable quantum systems . However, the analytical results suggest a simple scheme for numerical calculations and numerical simulations of the SRS process in terms of the six-vertex model of statistical mechanics .

The Hamiltonian of a quantum version of the standard SRS model can be written in the following form :

$$H=_{\mathrm{}}^{\mathrm{}}𝑑x\left\{i\overline{ϵ}_\sigma (x)\frac{}{x}ϵ_\sigma (x)+J\underset{a=1}{\overset{M}{}}\delta (xx_a)\overline{ϵ}_\sigma (x)\left[\sigma _{\sigma \mu }^+r_a^{}+\sigma _{\sigma \mu }^{}r_a^+\right]ϵ_\mu (x)\right\},$$
(1)
where the summation over repeated indexes is assumed. In Eq. (1), the operators of the “slow” envelopes of the laser (pump), $`E_L(x)`$, and Stokes, $`E_S(x)`$, components of light are combined into the isotopic spinor

$$ϵ_\sigma =\left(\begin{array}{c}E_L\\ E_S\end{array}\right),\overline{ϵ}_\sigma =(E_L^{},E_S^{})$$

with the commutator $`[ϵ_\sigma (x),\overline{ϵ}_\mu (y)]=\delta _{\sigma \mu }\delta (xy)`$, $`\sigma ,\mu =L,S`$. To keep our terminology simple, we use hereafter the term “spin” for the field’s subindex $`\sigma `$. Thus, the radiation field is described in terms of particles with spin up (a laser photon) and down (a Stokes photon).

To account for a possible saturation of a resonance transition of a medium, we treat the medium as a discrete set of $`M`$ two-level atoms with the coordinates $`x_a`$, $`a=1,\mathrm{},M`$, along the light propagation axis - the $`x`$ axis. All the atoms are assumed to be positioned on an interval of size $`l`$, $`0<x_1<\mathrm{}<x_M<l`$. As usual, they are described by the spin-$`\frac{1}{2}`$ operators $`𝐫_a=(r_a^x,r_a^y,r_a^z)`$, $`r^\pm =(r^x\pm ir^y)/2`$. The Pauli matrices $`\sigma _{\sigma \mu }^\pm `$ act in the spin space of the radiation. The first term of the model Hamiltonian describes a free propagation of the field components in the positive direction of the sample axis, while the second one represents a medium-field interaction with the coupling constant $`J`$, the two-level atoms playing the role of “spin impurities”.

In the classical limit and in the limit of a continuous description of the medium, the equations of motion for the dynamical variables of the field plus medium system coincide with the equations of the classical SRS theory . In the quantum approach, integrability of the SRS model is obvious, because the Hamiltonian (1) coincides with the Hamiltonian of the Kondo model with anisotropic exchange coupling whose integrability has been proven by Wiegmann .

The scattering of the $`j`$-th particle on the $`a`$-th impurity is described by the scattering matrix

$$S_{ja}=\frac{1}{2}(11+\sigma _j^zr_a^z)+\frac{b}{2}(11\sigma _j^zr_a^z)+c(\sigma _j^+r_a^{}+\sigma _j^{}r_a^+),$$
(2)
where $``$ implies the tensor multiplication, and the scattering amplitudes $`b=\mathrm{cos}J`$ and $`c=i\mathrm{sin}J`$ correspond to the elastic (Rayleigh) and inelastic (Stokes) channels, respectively. For what follows, it is convenient to use the independence of the $`S`$-matrix of a particle momentum  and to formulate the SRS problem only in terms of spin variables of particles and impurities.

The initial state of the spontaneous SRS problem is then given by

$$|\mathrm{In}=|\overline{\mathrm{\Omega }}_\text{p}|\mathrm{\Omega }_\text{i},$$
(3)
where the states of $`N`$ incident particles (all spins up) and $`M`$ impurities (all spins down) are defined as $`\sigma _j^+|\overline{\mathrm{\Omega }}_\text{p}=0`$ and $`r_a^{}|\mathrm{\Omega }_\text{i}=0`$. It is convenient to represent the particle-impurity scattering matrix (2) as a $`2\times 2`$ matrix acting in the spin space of a particle, while its matrix elements contain the spin variables of an impurity,

$$S=\left[\begin{array}{cc}\frac{1}{2}[(1+b)+(1b)r_a^z]& cr_a^{}\\ cr_a^+& \frac{1}{2}[(1+b)(1b)r_a^z]\end{array}\right].$$

The subsequent scattering of the $`j`$-th particle on the system of $`M`$ impurities is then described by the ordered scalar product of matrices $`S_{ja}`$,

$$L_j=S_{jM}\mathrm{}S_{j1}=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right],$$

where the operators $`A,B,C`$, and $`D`$ act in the spin space of impurities. Acting on the initial state (3), the monodromy matrix $`L_j`$ obviously reduces to the matrix

$$U_j=A+C\sigma _j^{},$$

where the first term corresponds to the elastic channel of photon scattering on the atomic system, while the second one represents the inelastic channel. In the latter, an incident laser photon is converted into a Stokes one, and simultaneously an excitation in the atomic system is created by the operator $`C`$. The exact final state of the spontaneous SRS problem is thus given by the ordered product of the matrices $`U_j`$,

$$|\mathrm{Out}_\mathrm{N}=(A+C\sigma _N^{})\mathrm{}(A+C\sigma _1^{})|\mathrm{In}.$$
(4)
Unlike the SRS problem in a pointlike geometry of the atomic system , the algebra of the elements of the monodromy matrix $`L`$ in an extended medium is not closed. Therefore, calculations of expectation values on the out-state (4) run into great mathematical troubles. Here, we compute the scattering amplitudes only for first two incident photons to clarify a microscopic quantum mechanism of the SRS phenomenon.

The first photon ($`j=1`$) propagates in the unexcited atomic system and the explicit expression for the out-state is easily found to be

$$|\mathrm{Out}_1=b^M|\mathrm{In}+c\underset{a=1}{\overset{M}{}}b^{(a1)}r_a^+|\mathrm{\Omega }_i\sigma _1^{}|\overline{\mathrm{\Omega }}_p.$$
(5)
Eq. (5) has a clear physical meaning. The first term corresponds to the elastic channel with the scattering amplitude $`𝒜_L^{(1)}=b^M`$ and the probability $`P_L^{(1)}=b^{2M}`$. The second term represents the inelastic channel, in which the laser photon is converted into a Stokes photon, and simultaneously an excitation in the atomic system is created. The probability amplitude of finding the $`a`$-th atom in the excited state, $`\varphi _a=cb^{(a1)}`$, is naturally treated as the wave function of a one-particle medium excitation created by the transmitted photon.

In this expression, the term $`c=i\mathrm{sin}J`$ is the probability amplitude to excite the $`a`$-th atom, while $`b^{(a1)}=(\mathrm{cos}J)^{(a1)}`$ is the probability amplitude that the scattering on the preceding $`a1`$ atoms was elastic. Since $`J1`$, the wave function of the excitation takes an obvious exponential form

$$\varphi (x)=iJ\mathrm{exp}(J^2\rho x/2).$$
(6)
Here, aa atomic number $`a`$ is represented as $`a=\rho x`$, where $`\rho =M/l`$ is the linear density of the number of atoms, and $`x`$ is the atomic coordinate measured from the left edge of the medium.

Let us consider now the scattering process of the second incident laser photon ($`j=2`$). This photon propagates in the medium whose quantum state has been prepared by the first photon. The probability of the elastic scattering of the second particle is obviously given by

$$P_L^{(2)}=\mathrm{Out}_1|AA|\mathrm{Out}_1,$$
(7)
and we have

$$A|\mathrm{Out}_1=b^Mb^M|\mathrm{In}+b^{(M1)}c\left[\underset{a=1}{\overset{M}{}}b^{(a1)}r_a^++c^2\underset{a_1=1}{\overset{M}{}}\underset{a_2=1}{\overset{a_11}{}}b^{(a_21)}r_{a_2}^+\right]|\mathrm{\Omega }_i\sigma _1^{}|\overline{\mathrm{\Omega }}_p.$$
(8)
Here, the first term describes the elastic propagation of the particle in the unexcited medium, while the second and third terms correspond to the particle propagation in the medium excited by the first incident photon. The second term represents a sub-channel of the elastic channel of scattering in which a particle does not change the direction of its spin propagating within the medium as a laser photon, $`1\stackrel{L}{}M`$. The third term represents the other possible sub-channel of the elastic scattering in which a particle changes the direction of its spin twice. Since one of the atoms, say $`a_1`$, is excited, a particle can change a spin direction first exciting an atom with a number $`a_2<a_1`$, and then deexciting the atom $`a_1`$, $`1\stackrel{L}{}a_2\stackrel{S}{}a_1\stackrel{L}{}M`$. Thus, on an interval between atoms $`a_2`$ and $`a_1`$ a particle propagates as a Stokes photon, but comes out of the medium as a laser photon.

Despite the factor $`c^2`$ being very small, the contribution of the third term to the probability $`P_L^{(2)}`$ is not small compared to the contribution of the second one. Moreover, omitting the third term, we immediately find that $`P_L^{(2)}>P_L^{(1)}`$, and hence the probability of the inelastic scattering falls with the growth of a photon number, $`P_S^{(2)}<P_S^{(1)}`$. The existence of two sub-channels of scattering for the second incident particle plays a crucial role for the nontrivial behavior of the probabilities of scattering. The interference of the scattering amplitudes in two possible sub-channels completely determines the growth of the Stokes component. In the first sub-channel, the excited state of the atomic system is not changed, while in the second sub-channel an excitation of the atomic system created by the first incident photon is shifted to the left edge of the medium.

One obtains from Eqs. (7) and (8),

$$P_L^{(2)}=p^{2M}+\left[4+p+\frac{M|c|^2(M|c|^24)}{1p^M}\right]p^{(M1)}(1p^M),$$
(9)
where $`p=b^2`$, and hence, in the limit $`J1`$,

$$\frac{P_L^{(2)}}{P_L^{(1)}}=12(MJ^2)^2+𝒪[(MJ^2)^4].$$
(10)
Thus, the interference of the scattering amplitudes in two sub-channels of the elastic channel of scattering drastically change the behavior of the probability $`P_L^{(2)}`$. Now, $`P_L^{(2)}<P_L^{(1)}`$, and hence, the probability of conversion of the second laser photon into the Stokes one, $`P_S^{(2)}=1P_L^{(2)}`$, is bigger than this probability for the first transmitted photon, $`P_S^{(1)}=1P_L^{(1)}`$, $`P_S^{(2)}>P_S^{(1)}`$.

The first derivative over time of the Stokes component intensity on the edge of an transmitted pulse, is easily found to be

$$\frac{d}{dt}I_S(t=0)=I_0\frac{P_S^{(2)}P_S^{(1)}}{P_S^{(1)}}2I_0(MJ^2)^2,$$
(11)
where $`I_0`$ is the intensity of an incident laser pulse. Thus, the intensity of the Stokes component at the initial moment of time grows with the velocity proportional to $`M^2`$; this is traditionally treated as a signature of the cooperative behavior of the SRS process.

Let us consider now an initial state of the scattering problem in which all incident particles are Stokes photons (all particle spins down), while all atoms are excited (all impurity spins up),

$$|\mathrm{In}=|\mathrm{\Omega }_\text{p}|\overline{\mathrm{\Omega }}_\text{i},$$
(12)
where as $`\sigma _j^{}|\mathrm{\Omega }_\text{p}=0`$ and $`r_a^+|\overline{\mathrm{\Omega }}_\text{i}=0`$. Then, the transmitted radiation contains obviously both the Stokes and laser components. It is clear also that the scattering process results in a deexcitation of atomic system.

Let us consider the case of a single excited atom. The probability amplitude $`A_{\text{ex}}(N)`$ to find the atom in the excited state after scattering $`N`$ incident laser photons is easily found from Eq. (2), $`A_{\text{ex}}(N)=(\mathrm{cos}J)^N`$. Taking into account that $`J1`$, we obtain for a time evolution of the probability amplitude

$$A_{\text{ex}}(t)=\left(1\frac{J^2}{2}\right)^N=\mathrm{exp}\left(\frac{1}{2}\gamma t\right).$$
(13)
Here, $`\gamma =J^2\rho `$, and the number of scattered photons at the moment of time $`t`$ is represented as $`N=\rho t`$, where $`\rho `$ is the linear density of incident photons.

As it should be expected, the excited atomic state decays exponentially with the decay constant $`\gamma `$ to the ground state. If

$$\rho \mathrm{},J0,\text{but}J^2\rho =\gamma =\text{constant},$$
(14)
as is already assumed in Eq. (13), an incident field is obviously eliminated from consideration, and the decay can be treated as spontaneous one, in which a single “Stokes” photon is emitted. In this limit, we can replace $`JE_L\sqrt{\gamma }`$ and $`E_S(x)E(x)`$ to reduce the Hamiltonian of the SRS problem (1) to the standard Hamiltonian of the quantum Maxwell-Bloch model,

$$H=i𝑑xE^{}(x)\frac{}{x}E(x)\sqrt{\gamma }\underset{a}{}𝑑x\left[E^{}(x)r_a^{}+r_a^+E(x)\right]\delta (xx_a).$$
(15)
An equivalence of these two models in the limit (14) has been used  to study the SRS problem in the framework of Maxwell-Bloch equation. It is clear that a solution of the more general SRS problem can be reduced to the SF problem only under some special conditions . However, a solution of the SF problem, including statistical properties of emitted photons , is exactly derived from a solution of the SRS problem in the limit (14).

I am thankful to S. John for stimulating discussions.