text:
# Anharmonicity, vibrational instability and Boson peak in glasses

## I Introduction

One of the most striking properties of glasses is a maximum in the inelastic scattering intensity observed in neutron and Raman scattering at frequencies between 0.5 and 2 THz, far below the Debye frequency . This so called Boson peak (BP) indicates an excess of low frequency vibrations over the Debye value which is given by the sound waves. It is seen as a maximum in the vibrational DOS divided by $`\omega ^2`$, $`G(\omega )/\omega ^2`$ but not necessarily in $`G(\omega )`$ itself . The BP is a general feature, its magnitude varies strongly between materials.

Despite numerous efforts, the BP remains one of the most intriguing problems of solid state physics. Some authors attribute the BP to vibrations of clusters of atoms of typical sizes . The physical origin of these clusters in homogeneous amorphous media remains unclear and they have not been identified in numerical simulations.

Another popular qualitative explanation of the Boson peak is a softening of acoustic phonons by static disorder  due to elastic Rayleigh scattering. However, even the most optimistic estimates show that the Rayleigh contribution is at least 4 times too small to explain the experimental data on thermal conductivity in glasses . This mechanism is also in contradiction to the linear dispersion law for acoustical phonons at the Boson peak frequency seen in molecular dynamics in Ref. .

Sometimes the BP is related to low lying optic modes of parental crystals . Whereas the BP is a general feature of glasses such crystal structures with soft optic modes cannot be identified always. Such a mechanism is possible in some cases. Yet it remains a puzzle how these crystalline peaks are transformed in the glassy state to a shoulder in the vibrational DOS.

Recent work on harmonic lattice models demonstrated that softening of disordered force constants can smear and push to low frequencies peaks which exist in the crystalline DOS . In another approach the vibrations of a random distribution of atoms, interacting with a Gaussian-shaped pair potential, was studied  in a harmonic scalar approximation. Reducing the density the system becomes unstable. Approaching this instability a low frequency peak appears in $`G(\omega )/\omega ^2`$ which resembles the BP. The main drawback of these models is their neglect of the static displacement of the atoms in response to disorder. In real glasses where short range order is conserved such feedback always occurs since the forces between the atoms strongly depend on their distance. In particular, the above models have no built in mechanism to stabilize vibrations with negative $`\omega ^2`$, unstable modes.

The proposed models of the Boson peak do not account for anharmonicity effects which, as we will show in the present paper, become very important especially for small force constants. Glasses, at low frequencies/temperatures, are highly anharmonic as seen in most of their macroscopic thermodynamic functions. Anharmonicity and static displacements, together, stabilize otherwise unstable vibrational modes. This does, however, not imply that the vibrations at the Boson peak are anharmonic. Anharmonicity is essential in forming the equilibrium structure, and thus the force constants which determine the proper harmonic spectrum of low frequency modes.

Another important point is that the previous explanations do not relate the Boson peak, which is one of the universal properties of glasses, to other universal properties, such as the two-level systems which dominate the low-temperature behavior or the plateau in the thermal conductivity at moderate temperatures. This relationship naturally emerges in our approach.

## II Quasi-Localized Vibrations (QLV)

We present a new universal mechanism for the formation of a BP in glasses, out of an originally flat DOS. This arises from the following 3 features: 1) quasi-local vibrations (resonant states) (QLV) with a smooth, structureless initial DOS, $`g_0(\omega )`$ at low $`\omega `$, 2) an elastic interaction between them and 3) stabilization by anharmonicity when the system becomes unstable because of interaction.

As the well known two level systems, the QLV are a typical feature of disordered systems. They are additional modes and are characterized by a large vibrational amplitude of some group of atoms. Their existence in glasses was predicted in Refs. . They can be described as low frequency harmonic oscillators (HO) which couple bilinearly to the sound waves, see Refs. . This in turn leads to a dipole-dipole interaction between different HO. The importance of the elastic interaction between local defects in glasses has been stressed by Yu and Leggett  and Grannan et al. .

The microscopic origin of QLV in disordered systems varies and depends on the type of disorder. Broadly speaking, the QLV can be divided into two groups. First there are materials where QLV exist independently of the structural disorder typical for structural glasses and amorphous solids. Orientationally disordered (“plastic”) crystals belong to this group . In these materials some molecular groups librate with low frequencies. In harmonic approximation these soft librations can be identified with QLV. The local potentials for the librational motion can vary from site to site. In this case we have a distribution of the librational frequencies. The librations couple to the sound waves which in turn induces an interaction between them . It depends on details of the material whether the interaction is strong enough to reconstruct the original spectrum of QLV completely or only partially. This effect is seen in recent measurements of dielectrtric loss spectra of ortho-carborane . Similarly coordination defects in covalent materials can lead to QLV as was observed by Biswas et al.  in a simulation of amorphous Si (see also recent works  and ). QLV can originate from numerous defects such as off center ions or interstitial atoms. Depending on the “size” of the defect the QLV involve more or fewer atoms. Interstitial atoms are the prototype of a topological point defect. QLV of interstitial atoms in FCC metals were studied extensively in the past, see Ref.  for a review. These QLV have effective masses of four atomic masses and the crystal structure is strongly distorted by the defect and the low frequency of the librational QLV can be traced to the local strain. Already low concentrations of these interstitials are sufficient to destroy the crystalline structure completely. This is utilized in the interstitialcy model of glass formation .

This leads directly to the second group where the QLV result directly from disorder. Such modes are regularly found in computer simulations, e.g. for soft spheres , SiO2 , Se , Ni-Zr , Pd-Si and Au-Si , NiB , in amorphous ice , in amorphous and quasicrystalline Al-Zn-Mg  and in simple dense fluids .

In these simulations the QLV were observed as localized vibrations with frequencies below the minimal sound wave frequency allowed by the size of the simulated sample. A simple indicator of these “size-localized” QLV is the scaling of the participation ratio with system size. Increasing the system size the minimal sound wave frequency drops and the QLV are no longer localized in the simulation but show the typical properties of resonant modes, i.e. of low frequency local vibrations which couple bilinearly to the sound waves. The exact eigenvectors of the interacting system of QLV and sound waves are superpositions of these two types of modes. The effect of system size on the appearance of QLV in simulations was discussed in detail in Ref. . There it was shown that the exact eigenvectors at frequencies up to and above the boson peak can be decomposed into extended sound wave like modes and the local cores of QLV. The latter correspond to the harmonic oscillators of this work.

The physical origin of these disorder induced QLV can be traced to local irregularities of the amorphous structure. In dense packed metallic glasses these originate e.g. from the conflict of the local dense packing (icosahedral packing) and global dense packing (fcc or hcp) . These local irregularities can be seen as centers of local strains .

We expect such local strain centers to be ubiquitous. The strains will have broad distributions, which will lead to broad distributions of QLV frequencies. Whereas local strains and QLV will be a general property of glasses the atomistic structure of QLV reflects the structure of the considered material. In dense packed metallic systems the cores of the QLV have been found to be chain-like . In SiO2 they are formed by a coupled rotations of SiO4 tetrahedra . In Se one has coupled chains and rings , etc.

Another possible (and natural) mechanism of the QLV formation are low lying optical modes in parental crystals. Disorder in amorphous material would destroy the long range coherence of optical modes. This makes them practically indistinguishable from quasi-local modes.

Together with the tunneling systems, the QLV form the main ingredient of the soft potential model  (see review ). They manifest themselves in experimental values, e. g. the excess specific heat  and the plateau in the thermal conductivity , in inelastic light  and neutron scattering  and they are observed in many numerical calculations, cited above.

## III Density of states of Quasi-Localized Vibrations (QLV)

One might think a priori that the QLV can have an arbitrary DOS, $`g(\omega )`$, depending on the particulars of the glass. We will show, that due to the interaction between the HO, $`g(\omega )`$ is a universal function at low frequencies. This universality stems from the vibrational instability of the spectrum which occurs in nearly all systems of interacting HO. Anharmonicity stabilizes the system in new minima, and thus reconstructs the DOS to a new harmonic spectrum. As a result $`g(\omega )/\omega ^2`$ acquires a maximum without a peak in $`g(\omega )`$ itself.

Below we will derive the following form of the reduced DOS of these harmonic resonant modes (excluding the Debye part at low $`\omega `$ which is not seen in Raman scattering)

$`{\displaystyle \frac{g(\omega )}{\omega ^2}}={\displaystyle \frac{3C}{\pi \omega ^{}}}\left({\displaystyle \frac{\omega ^{}}{\omega }}\right)^4\left[z_1^2(\omega )+z_2^2(\omega )\right]^1`$ (1)
$`\times \left[{\displaystyle \frac{1}{2z_1(\omega )}}\mathrm{ln}{\displaystyle \frac{z_1(\omega )+1}{z_1(\omega )1}}+{\displaystyle \frac{1}{z_2(\omega )}}\mathrm{tan}^1{\displaystyle \frac{1}{z_2(\omega )}}\right]`$ (2)
where $`C`$ is a constant (Eq. (23)) and

$$z_{1,2}(\omega )=\frac{1}{2}\sqrt{\sqrt{9+8(\omega ^{}/\omega )^6}\pm 3}.$$
(3)
The function $`g(\omega )/\omega ^2`$ is plotted in Fig. 1. It depends on a single parameter, $`\omega ^{}`$ characterizing the position of the Boson peak. The maximum of $`g(\omega )/\omega ^2`$, the Boson peak, is at $`\omega _b1.1\omega ^{}`$. For small frequencies, $`\omega \omega _b`$, $`g(\omega )\omega ^4`$ while for large ones, $`\omega \omega _b`$, $`g(\omega )\omega `$.

## IV Vibrational instability

To illustrate our central idea of a vibrational instability, we start with a pair of interacting HO immersed in an elastic continuum. The potential energy is given by

$$U_{\mathrm{har}}(x_1,x_2)=M_1\omega _1^2x_1^2/2+M_2\omega _2^2x_2^2/2I_{12}x_1x_2.$$
(4)
Here $`x_{1,2}`$ are the HO coordinates, $`M_{1,2}`$ the masses and $`\omega _{1,2}`$ the bare frequencies of the two HO, i. e. neglecting the bilinear interaction. The interaction strength is given by

$$I_{12}=g_{12}J/r_{12}^3,J\mathrm{\Lambda }^2/\rho v^2$$
(5)
where $`g_{12}`$ accounts for the relative orientation of the HO, $`r_{12}`$ is their distance, $`\rho `$ is the mass density of the glass and $`v`$ is a sound velocity. The interaction between the HO is due to the coupling between a single HO and the surrounding elastic medium (the glass). This HO-phonon coupling has the form

$$_{\mathrm{int}}=\mathrm{\Lambda }x\epsilon ,$$

where $`\mathrm{\Lambda }`$ is the coupling constant and $`\epsilon `$ the strain.

Diagonalization of Eq. (4) yields two frequencies

$$\stackrel{~}{\omega }_{1,2}^2=\frac{\omega _1^2+\omega _2^2}{2}\sqrt{\left(\frac{\omega _1^2\omega _2^2}{2}\right)^2+\frac{I_{12}^2}{M_1M_2}}.$$
(6)
The smaller value, $`\stackrel{~}{\omega }_1^2`$, becomes negative when the interaction $`I_{12}I`$ exceeds the threshold (critical) value

$$I_c\omega _1\omega _2\sqrt{M_1M_2}.$$
(7)
This instability persists also in a system of many interacting HO. In a real physical system, anharmonic forces always stabilize an embedded HO in a nearby minimum of the potential energy. The position of this minimum depends on the interaction between HO. We are thus confronted with the many-body problem of finding the minima of the potential energy for a system of interacting anharmonic oscillators, similar to the one considered in Ref. . The new frequencies, in these new minima, are real and different from the original ones. The harmonic vibrational spectrum is reconstructed. We will call this anharmonicity limited vibrational instability.

## V Stabilization by anharmonicity

We will now show that for weak interaction, $`I`$, the reconstructed DOS has, below a characteristic frequency $`\omega _c|I|`$,a universal form irrespective of its original form. First, due to interaction, it becomes a linear function of frequency, $`g(\omega )\omega `$. Secondly, the displacements of the previously unstable oscillators from their old equilibrium positions create static random forces which cause a second reconstruction of the DOS below another frequency $`\omega _b\omega _c`$. Due to so called sea-gull singularity  at $`\omega =0`$ the linear DOS is reconstructed to $`g(\omega )\omega ^4`$ for $`\omega <<\omega _b`$. Together, these two reconstructions produce a maximum of $`g(\omega )/\omega ^2`$ at $`\omega =\omega _b`$.

Let us consider a number of randomly distributed, interacting HO with concentration $`n_0`$ and an initial DOS, $`g_0(\omega )`$ (normalized to unity), in the frequency range from 0 to $`\omega _0`$, where $`g_0(\omega )`$ is a monotonously increasing function of $`\omega `$. For the harmonic part of the interaction we take the generalization of Eq. (4) and add an anharmonic term to stabilize the system

$$U_{\mathrm{anhar}}=(1/4)\underset{i}{}A_ix_i^4,A_i>0.$$
(8)
We will take the interaction $`I`$ to be the small parameter of our theory, i.e. we assume that the typical random interaction, $`I`$, between neighboring HO is much smaller than the typical values of $`M\omega _0^2`$. As $`|I|M\omega _0^2`$ frequencies of order $`\omega _0`$ will be practically unaffected by the interaction whereas HO with frequencies $`\omega <\omega _c`$ will be displaced to new minima, where

$$\omega _c|I|/M\omega _0\omega _0.$$
(9)
Since the concentration of unstable HO is much smaller than the one of the stable ones a low frequency oscillator is typically surrounded by high frequency ones. We can simplify our consideration by again considering pairs of HO, one with a low frequency $`\omega _1\omega _c`$ and the other one from the bulk of the HO with a frequency $`\omega _2`$ of the order of $`\omega _0`$ (see Appendix A for the general case). Due to the combined action of interaction and anharmonicity the two HO will be displaced into new minima, $`x_{10}`$ and $`x_{20}`$, given by the equations ($`I_{12}I`$)

$$\begin{array}{ccc}\hfill Ix_{20}& =& x_{10}(M_1\omega _1^2+A_1x_{10}^2),\hfill \\ \hfill Ix_{10}& =& x_{20}(M_2\omega _2^2+A_2x_{20}^2).\hfill \end{array}$$
(10)
For $`|I|>I_c`$ we need the nonzero solutions of these equations. Expanding around either minimum we find the new (harmonic) frequencies from the secular equation

$$\left|\begin{array}{cc}\alpha _1M_1\omega ^2& I\\ I& \alpha _2M_2\omega ^2\end{array}\right|=0$$
(11)
with

$$\alpha _i=M_i\omega _i^2+3A_ix_{i0}^2,i=1,2.$$
(12)
¿From the condition $`\omega _1\omega _2`$ follows $`x_{20}x_{10}`$ and, therefore, the term $`A_2x_{20}^2`$ in Eq. (10) can be neglected giving

$$x_{20}=(I/M_2\omega _2^2)x_{10}.$$
(13)
and

$$x_{10}=\omega _1\sqrt{M_1/A_1}\sqrt{(I/I_c)^21}.$$
(14)
As a result we get from (11) under the condition $`\omega _1/\omega _21`$ with this accuracy the new frequencies $`\stackrel{~}{\omega }_2=\omega _2`$ and

$$\stackrel{~}{\omega }_1^2=2\omega _1^2\left[(I/I_c)^21\right].$$
(15)
The smaller frequency (15) is the solution of the linear equation (compare with Eq. (11))

$$M_2\omega _2^2(\alpha _1M_1\omega ^2)=I^2.$$
(16)
It is remarkable that for weak interaction the strength of the anharmonicity $`A_i`$ does not enter the renormalized frequency (15).

Near the threshold where $`(|I|I_c)/I_c1`$, the smaller frequency squared $`\stackrel{~}{\omega }_1^2`$ is proportional to $`(|I|I_c)/I_c`$. Provided the distribution of the random quantity $`I`$ is smooth one gets, therefore, below $`\omega _c`$ a linear DOS ($`\stackrel{~}{g}(\omega )\omega `$) irrespective to the initial form of $`g_0(\omega )`$. In Appendix A it is shown that the same result holds if one has a low frequency HO surrounded by several high frequency ones. Our numerical calculations (see Section VII) also show that this case is typical.

## VI The Boson peak

If the low-frequency HO with their reconstructed linear DOS were isolated, the problem would be solved. There is, however, a further interaction between these oscillators which we have not taken into account so far. The low-frequency HO, displaced from their equilibrium positions, create random static forces $`f`$. The force $`f_i`$ exerted on the $`i`$th oscillator by the $`j`$th one is

$$f_i=I_{ij}x_{j0}.$$
(17)
In a purely harmonic case, these linear forces would not affect the frequencies. Anharmonicity, however, renormalizes the low frequency part of the spectrum , a manifestation of the sea-gull singularity. Consider an anharmonic oscillator under the action of a random static force $`f`$

$$U(x)=Ax^4/4+M\omega _1^2x^2/2fx$$
(18)
where $`\omega _1`$ is the oscillator frequency in the harmonic approximation. The force $`f`$ shifts the equilibrium position from $`x=0`$ to $`x_00`$, given by

$$Ax_0^3+M\omega _1^2x_0f=0,$$
(19)
where the oscillator has a new (harmonic) frequency

$$\omega _{\mathrm{new}}^2=\omega _1^2+3Ax_0^2/M.$$
(20)
If $`\stackrel{~}{g}_1(\omega _1)`$ is the distribution function of frequencies $`\omega _1`$ and $`P(f)`$ is the distribution of random forces, then the renormalized DOS is given by

$$g(\omega )=\underset{0}{\overset{\mathrm{}}{}}\stackrel{~}{g}_1(\omega _1)𝑑\omega _1\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑fP(f)\delta \left(\omega \omega _{\mathrm{new}}\right).$$
(21)
As the forces between the HO are proportional to $`r_{ij}^3`$ their sum is Lorentzian distributed (see Appendix B):

$$P(f)=\frac{1}{\pi }\frac{\delta f}{f^2+(\delta f)^2}.$$
(22)
Assuming $`\omega \omega _c`$ and integrating Eq. (21) with $`\stackrel{~}{g}_1(\omega _1)=C\omega _1`$ we arrive at the integral

$$\frac{g(\omega )}{\omega ^2}=\frac{6C}{\pi \omega ^{}}\left(\frac{\omega }{\omega ^{}}\right)^2\underset{0}{\overset{1}{}}\frac{dt}{1+\left(\omega /\omega ^{}\right)^6t^2(32t^2)}$$
(23)
with

$$\omega ^{}=\sqrt{3}A^{1/6}(\delta f)^{1/3}/\sqrt{M}$$
(24)
and after integration finally Eq. (1) is obtained.

For small frequencies, below the Boson peak, $`\omega \omega _b`$ only small random forces $`f`$ contribute to the second integral in Eq. (21). In this case the distribution function $`P(f)`$ can be approximated by a constant value, $`P(0)`$, and we get from (21)

$$g(\omega )\omega ^3\underset{0}{\overset{\omega }{}}𝑑\omega _1\frac{\omega _1}{\sqrt{\omega ^2\omega _1^2}}\omega ^4.$$
(25)
As a result at low frequencies the renormalized DOS, $`g(\omega )\omega ^4`$, Ref. . For sufficiently large frequencies, $`\omega \omega _b`$ the action of random static forces on the HO spectrum can be discarded. In this case we recover the linear DOS, $`g(\omega )\omega `$.

The BP frequency, $`\omega _b1.1\omega ^{}`$, is determined by the characteristic value of the random static force $`\delta f`$, Eq. (17), acting on an HO with the characteristic frequency $`\omega _c`$. According to Eq. (17), it is due to the interaction between HO with frequencies of order of $`\omega _c`$, i.e.

$$I_{ij}^{(c)}Jn_c,JI/n_0M\omega _c/n_0g_0(\omega _0)$$

where $`n_cn_0g_0(\omega _c)\omega _c`$ is the concentration of these HO. The characteristic displacement of a low-frequency HO (14) from the equilibrium position is $`x_{j0}\omega _c\sqrt{M/A}`$. As a result, we get the estimate

$$\delta fM\sqrt{\frac{M}{A}}\omega _c^3\frac{g_0(\omega _c)}{g_0(\omega _0)}$$
(26)
and according to (24)

$$\omega _b\omega _c\left[g_0(\omega _c)/g_0(\omega _0)\right]^{1/3},\omega _b\omega _c.$$
(27)
Again, in lowest order the anharmonicity $`A`$ does not enter this formula.

As a result we get a following estimate for the reconstructed DOS

$$g(\omega )\{\begin{array}{cc}g_0(\omega ),\hfill & \omega >\omega _c,\hfill \\ \omega g_0(\omega _c)/\omega _c,\hfill & \omega _b<\omega <\omega _c,\hfill \\ \omega ^4g_0(\omega _0)/\omega _c^4,\hfill & \omega <\omega _b.\hfill \end{array}$$
(28)
If the DOS of the noninteracting oscillators is given by a power law, $`g_0(\omega )\omega ^n`$, the BP frequency $`\omega _b`$ scales with the interaction strength $`I`$ as

$$\omega _b|I|^{1+n/3}.$$
(29)
Since in accordance with (9) $`\omega _c|I|`$, we have always $`\omega _b\omega _c`$ for $`n>0`$.

## VII Numerical simulation

To test our ideas by numerical simulations, we placed $`N`$ oscillators with frequencies $`0<\omega _i<1`$ on a simple cubic lattice with lattice constant $`a=1`$ and periodic boundary conditions. To simulate random orientations of the oscillators we took for $`g_{ij}`$, Eq. (5), random numbers in the interval $`[0.5,0.5]`$. The masses, $`M_i`$, and anharmonicity parameters, $`A_i`$, were put to 1. The DOS for the noninteracting oscillators was taken as $`g_0(\omega )\omega ^n`$, with $`n=1,2,3`$.

Using the potential energy given by the generalization of Eq. (4) plus the anharmonicity (Eq. (8)) we then minimized the potential energy, and in the usual harmonic expansion around this minimum calculated the DOS for different interaction strengths, $`J`$. This was repeated for up to 10000 representations. To check for size dependence we did the calculations for different $`N`$. Apart from the case $`J=0.07`$ the results did not change between $`N=2097`$ and $`N=4096`$.

The predicted change-over in the $`\omega `$-dependence of the DOS at two characteristic frequencies $`\omega _c`$ and $`\omega _b`$ and the linear part in between can be clearly observed in a log-log representation, Fig. 2, for $`g_0(\omega )\omega ^2`$ and $`J=0.1`$. ¿From the calculated eigenvectors we find that, as expected, at the lowest frequencies the HO are weakly coupled whereas near and above $`\omega _b`$ the eigenmodes are complicated superpositions of many HO.

Fig. 3 shows the dependence of the simulated $`g(\omega )/\omega ^2`$ on the interaction strength $`J`$. We can see the general increase of $`\omega _b`$ and related decrease of the BP intensity with increasing $`J`$. Our simulations cover one decade in BP frequencies. The insert shows that, in full agreement with our predictions (see Eq. (9) and Eq. (29)), the crossover frequencies change with interaction $`J`$ as $`\omega _cJ`$ and $`\omega _bJ^{1+n/3}`$.

## VIII Discussion and comparison with experiment

In Fig. 4 we compare our theoretical curve, Eq. (1), with Raman scattering data of lithium borate glasses  with different compositions. The agreement is remarkably good over the whole composition range. This supports the idea of a universal shape of the Boson peak . The shift of the BP to higher frequencies with increasing concentration of Li2O can be explained by an increase of the total concentration of QLV and consequently of their interaction.

The agreement between theory and experiment is not confined to this class of material. This is exemplified by Fig. 5 which shows a comparison of the theoretical curve with numerical simulation results and neutron and Raman scattering data for 3 different glasses.

One of the most important results of our theory is the predicted linear frequency dependence of the density of vibrational states above the Boson peak. It stems from the vibrational instability of interacting harmonic modes. Such linear behavior has been observed in many numerical simulations on different glasses and model disordered systems . It is also in a good agreement with many experimental results  where the vibrational DOS has a section with near linear frequency dependence. Fig. 6 shows this for vitreous silica at different temperatures. Above the Boson peak the DOS increases approximately linearly with frequency.

In this paper we dealt with the case of weak interaction between HO. If the interaction is increased the characteristic frequencies $`\omega _b`$ and $`\omega _c`$ grow and the gap between them narrows and finally disappears. Then our ”BP” in $`g(\omega )/\omega ^2`$ superimposes the “boundary peak” in $`g_0(\omega )`$ at the edge of the assumed spectrum of naked (noninteracting) QLV. The BP can no longer be distinguished from the boundary peak or from a possible equivalent maximum in $`g_0(\omega )`$. It is possible that in some cases, e.g. orientational glasses, $`g_0(\omega )`$ has a pronounced peak which is still visible after the reshaping of the DOS by interaction.

Similar models with strong coupling between oscillators were previously investigated by molecular dynamic simulations  and by the replica method . For example in the Ref.  Grannan et al. assumed that the dynamics ”is completely dominated by the interaction between defects mediated by the strain field”. The local field (harmonic and anharmonic contribution) was neglected. We consider the opposite case where the local field (harmonic and anharmonic) dominates and interaction between dipoles is weak.

## IX Conclusion

In conclusion we present a universal picture of the BP formation in glasses. We have shown that the low frequency quasi-localized harmonic modes in glasses are destabilized by the weak bilinear interaction between them. Anharmonicity stabilizes the system in a new minimum of configuration space. It completely reconstructs the low frequency part of the spectrum (at $`\omega <\omega _c`$) and the Boson peak feature (at $`\omega _b\omega _c`$) naturally emerges. The thus created boson peak has a material independent shape. At low frequencies, below the BP, the vibrational DOS increases as $`g(\omega )\omega ^4`$, and above the BP, as $`g(\omega )\omega `$.

Although the anharmonicity is responsible for this effect, the final spectrum of stable vibrations remains harmonic. A remarkable feature of the presented theory is that the strength of the anharmonicity does not enter the new stable spectrum at all. It looks as if the anharmonicity does all the work, it stabilizes the system in a new minimum and reconstructs the spectrum and then, like a ghost, disappears. Therefore, the discussed phenomenon is independent of the variation of the anharmonicity between different materials. The only parameter entering the final density of states is the strength of interaction between the HO.

Compared to previous work, a main new result of our approach is the natural emergence of the BP on the unstructured, flat low frequency part of the initial spectrum $`g_0(\omega )/\omega ^2`$ where the DOS previously had no peaks. For small interactions the BP frequency is much smaller than the Debye frequency value. It shifts with interaction strength, $`I`$, which explains the large variety of BP magnitudes found in experiment. In contrast with previous models, a natural connection of the Boson peak phenomena with other universal properties of glasses is established.

The authors VLG and DAP gratefully acknowledge the financial support of the German Ministry of Science and Technology and the hospitality of the Forschungszentrum Jülich where part of the work was done.

## A Cluster approach

Consider a cluster containing a low frequency oscillator with frequency $`\omega _1\omega _c`$ surrounded by a large number, $`s1`$, of HO with much higher frequencies $`\omega _j\omega _0`$. Inclusive their interaction, the total potential energy of the cluster is

$$U_{\mathrm{tot}}=\underset{i}{}\frac{M_i\omega _i^2}{2}+\frac{1}{4}\underset{i}{}A_ix_i^4\frac{1}{2}\underset{i,ji}{}I_{ij}x_ix_j.$$
(A1)
The equilibrium positions of the HO, $`x_{i0}`$, are given by the system of $`s`$ nonlinear equations

$$M_i\omega _i^2x_{i0}+A_ix_{i0}^3=\underset{ji}{}I_{ij}x_{j0},i=1,2,\mathrm{},s.$$
(A2)
In the case of instability ($`x_{i0}0`$), in analogy to the previously considered case of a pair of oscillators, the static displacements of the high-frequency oscillators are much smaller than the one of the low-frequency oscillator, $`x_{10}`$. Therefore, in leading order

$$x_{i0}=(I_{1i}/M_i\omega _i^2)x_{10},i1.$$
(A3)
Inserting these values into Eq. (A2) for $`i=1`$ we get

$$M_1\omega _1^2x_{10}+A_1x_{10}^3=x_{10}\underset{i1}{}\frac{I_{1i}^2}{M_i\omega _i^2}.$$
(A4)
Under the condition

$$M_1\omega _1^2<k\text{where}k\underset{i1}{}\frac{I_{1i}^2}{M_i\omega _i^2}$$
(A5)
the cluster becomes unstable and the low frequency oscillator is displaced to a new minimum

$$x_{10}=\sqrt{(kM_1\omega _1^2)/A_1}.$$
(A6)
In the opposite case, $`M_1\omega _1^2>k`$, the cluster is stable and $`x_{i0}=0`$.

The eigenfrequencies of the interacting oscillators are the solutions of the secular equation of order of $`s`$

$$\left|\begin{array}{cccc}\alpha _1M_1\omega ^2& I_{12}& \mathrm{}& I_{1s}\\ I_{21}& \alpha _2M_2\omega ^2& \mathrm{}& I_{2s}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ I_{s1}& I_{s2}& \mathrm{}& \alpha _sM_s\omega ^2\end{array}\right|=0.$$
(A7)
Here the $`\alpha _i`$ are given by Eq. (12). In leading order in $`I_{ij}/M\omega _0^2`$ the secular equation is a linear equation for $`\omega ^2`$ (compare Eq. (16))

$$(\alpha _1M_1\omega ^2)\underset{j1}{}M_j\omega _j^2\underset{i1}{}I_{1i}^2\underset{j1,i}{}M_j\omega _j^2=0$$
(A8)
or

$$\alpha _1M_1\omega ^2=\underset{i1}{}\frac{I_{1i}^2}{M_i\omega _i^2}=k$$
(A9)
and the new low frequency of the system of coupled oscillators is given by

$$\stackrel{~}{\omega _1}^2=\{\genfrac{}{}{0pt}{}{{\displaystyle \frac{1}{M_1}}(M_1\omega _1^2k),k<M_1\omega _1^2}{{\displaystyle \frac{2}{M_1}}(kM_1\omega _1^2),k>M_1\omega _1^2.}$$
(A10)
As in the case of a pair of oscillators, the anharmonicity has been used in the derivation of Eq. (A10) but does not appear in this or our final result, Eq. (1).

To derive of the new (reconstructed) DOS the distribution of $`k`$, $`\rho (k)`$, has to be calculated. Inserting Eqs. (5) and (A5) into the definition of $`\rho (k)`$ gives

$$\rho (k)=\delta \left(k\frac{J^2}{M}\underset{j1}{}\frac{g_{1j}^2}{r_{1j}^6\omega _j^2}\right).$$
(A11)
Here the angular brackets denote averaging over the positions of the $`s1`$ high frequency HO, their frequencies and orientations. For simplicity we take equal masses $`M_j=M`$ and for $`g_{ij}`$ a uniform distribution in the interval \[-1/2, 1/2\].

Using the Holtsmark method  (see Appendix B) one gets

$$\rho (k)=\frac{1}{\sqrt{2\pi }}\frac{B}{k^{3/2}}\mathrm{exp}\left(\frac{B^2}{2k}\right)$$
(A12)
where

$$B=\frac{\pi }{3}\sqrt{\frac{\pi }{2}}\frac{Jn_0}{\sqrt{M}}\frac{1}{\omega }_0\omega _c\sqrt{M}.$$
(A13)
Here $`n_0`$ is the total concentration of HO in the cluster and $`1/\omega _0`$ is the $`\omega ^1`$ moment of the normalized initial DOS $`g_0(\omega )`$. This formula is a more accurate definition of the characteristic frequency $`\omega _c`$ introduced in Eq. (9). Note that the distribution $`\rho (k)`$ (A12) belongs to an important class of one-sided stable distributions, Ref. .

Due to the combined action of interaction and anharmonicity the DOS is reconstructed to $`\stackrel{~}{g}(\omega )=2\omega \stackrel{~}{G}(\omega ^2)`$ with

$`\stackrel{~}{G}(\omega ^2)=\delta \left(\omega ^2\stackrel{~}{\omega }_1^2\right)_{k,\omega _1}`$ (A14)
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k\rho (k){\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑\omega _1^2G_0(\omega _1^2)\delta \left(\omega ^2\stackrel{~}{\omega }_1^2\right)`$ (A15)
and $`G_0(\omega _1^2)g_0(\omega _1)/2\omega _1`$. Using Eq. (A10) and integrating Eq. (A15) we obtain

$`\stackrel{~}{G}(\omega ^2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k\rho \left(k+{\displaystyle \frac{1}{2}}M\omega ^2\right)G_0\left({\displaystyle \frac{k}{M}}\right)+`$ (A17)
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑k\rho (k)G_0\left(\omega ^2+{\displaystyle \frac{k}{M}}\right).`$
For low frequencies, $`\omega \omega _c`$, $`\stackrel{~}{G}(\omega ^2)=\text{const}`$ and $`\stackrel{~}{g}(\omega )\omega `$, i.e. the reconstructed DOS is a linear function of $`\omega `$. For high frequencies the first term in Eq. (A17) can be discarded and the original DOS is reproduced, $`\stackrel{~}{G}(\omega ^2)=G_0(\omega ^2)`$ for $`\omega \omega _c`$.

## B Holtsmark method

### 1 Distribution of random forces

Let $`x_i`$ be a random value with zero mean, $`x=0`$, and finite $`|x|`$ and let $`𝐫_i`$ ($`N\mathrm{}`$) be Poisson-distributed random points in 3-d space with concentration $`n_c`$. The distribution function, $`P(f)`$, of the random values

$$f=\underset{i}{}\frac{x_i}{r_i^3}$$
(B1)
can then be calculated by the Holtsmark method .

$$P(f)=\delta \left(f\underset{i}{}\frac{x_i}{r_i^3}\right)\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\tau e^{if\tau }F(\tau )$$
(B2)
with

$$F(\tau )=\mathrm{exp}\left[i\tau \underset{i}{}\frac{x_i}{r_i^3}\right],$$
(B3)
where angular brackets denote averaging over $`x_i`$ and $`𝐫_i`$. Since the values $`x_i/r_i^3`$ are independent of each other

$$F(\tau )=e^{i\tau x/r^3}^N=\left(1\frac{1}{V}d^3𝐫1e^{i\tau x/r^3}_x\right)^N$$
(B4)
which in the limit $`N\mathrm{}`$ becomes

$$F(\tau )=\mathrm{exp}\left[n_cd^3𝐫1e^{i\tau x/r^3}_x\right].$$
(B5)
Using $`x=0`$ and changing the integration variable to $`y=|\tau ||x|/r^3`$ we get

$$F(\tau )=\mathrm{exp}\left[\frac{4\pi n_c}{3}|\tau ||x|\underset{0}{\overset{\mathrm{}}{}}\frac{dy}{y^2}(1\mathrm{cos}y)\right].$$
(B6)
The integral equals $`\pi /2`$ and Eq. (B2) is the Fourier transform of a Lorentzian distribution of random forces

$$P(f)=\frac{1}{\pi }\frac{\delta f}{f^2+(\delta f)^2}$$
(B7)
where the width of the distribution is given by

$$\delta f=\frac{2\pi ^2n_c}{3}|x|.$$
(B8)
### 2 $`\rho (k)`$ distribution

The same method can be applied to calculate the distribution of the random quantity $`k`$

$$k=\underset{i}{}\frac{x_i^2}{r_i^6\omega _i^2}.$$
(B9)
where the $`x_i`$ are random and uniformly distributed, $`x_0/2<x_i<x_0/2`$, $`𝐫_i`$ are ($`N\mathrm{}`$) Poisson-distributed random points in 3-d space (concentration $`n_0`$) and $`\omega _i`$ are random frequencies of HO distributed in the interval

$$0<\omega _i<\omega _0$$
(B10)
with a DOS, $`g_0(\omega )`$, normalized to unity.

Analogously to (B2) the distribution $`\rho (k)`$ can be written as

$$\rho (k)=\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\tau e^{ik\tau }K(\tau )$$
(B11)
with

$$K(\tau )=\mathrm{exp}\left[i\tau \underset{i}{}\frac{x_i^2}{r_i^6\omega _i^2}\right].$$
(B12)
Following the steps of the previous subsection we can write

$$K(\tau )=\mathrm{exp}\left[n_0d^3𝐫1e^{i\tau x^2/r^6\omega ^2}_{x,\omega }\right]$$
(B13)
and, introducing the new variable $`y=(x^2|\tau |)/(r^6\omega ^2)`$,

$$K(\tau )=\mathrm{exp}\left[\frac{2\pi }{3}n_0\frac{|x|}{\omega }_{x,\omega }\sqrt{|\tau |}(\alpha +i\beta \mathrm{sign}\tau )\right]$$
(B14)
where

$$\alpha =\underset{0}{\overset{\mathrm{}}{}}\frac{dy}{y^{3/2}}(1\mathrm{cos}y),\beta =\underset{0}{\overset{\mathrm{}}{}}\frac{dy}{y^{3/2}}\mathrm{sin}y$$
(B15)
It is straightforward to show that $`\alpha =\beta =\sqrt{2\pi }`$. Therefore

$$K(\tau )=\mathrm{exp}\left[B\sqrt{|\tau |}(1+i\mathrm{sign}\tau )\right]$$
(B16)
with

$$B=\frac{\pi }{3}\sqrt{\frac{\pi }{2}}n_0x_0\frac{1}{\omega }$$
(B17)
After integration in (B11) with $`K(\tau )`$ from (B16) we finally get

$$\rho (k)=\frac{1}{\sqrt{2\pi }}\frac{B}{k^{3/2}}\mathrm{exp}\left(\frac{B^2}{2k}\right).$$
(B18)


format: text

text:
# Rate-equation approach to atomic-laser light statistics.

## I Introduction

Interest in the statistics of light emitted by atomic lasers has been recently revived as a result of the realization of micro-lasers An et al. (1994). The main purpose of this paper is to show that expressions derived from rate equations coincide with those previously derived from Quantum Optics. This is so even when the generated light exhibits sub-Poissonian statistics. An introduction to that method as it pertains to sub-Poissonian light generation can be found in tutorial papers Arnaud (1995, 2001). Let us recall that light is called sub-Poissonian when the variance of the number of photo-detection events over some large time duration is less than the average number of events. Equivalently, one may say that the spectral density of the photo-current is less than the shot-noise level at small Fourier (or baseband) frequencies. Analytical expressions are obtained from rate equations in a straightforward manner. Namely, the expression for the internal cavity statistics of many 4-level atoms with a negligible spontaneous decay previously given in Eq. (4) by Ritsch et al. (1991) is recovered exactly (see Eq. (20) of the present paper). Similarly, 3-levels atoms expressions obtained by Khazanov et al. (1990) are recovered. But coherently-pumped 3-level atoms lasers were apparently not treated earlier. If the upper and lower decay times of 4-level atoms tend to zero, the laser is equivalent to a 2-levels atoms laser with Poissonian pump Chusseau et al. . In that limit, the expressions reported in Arnaud and Estéban (1990) in 1990 are recovered. When the above approximations are not applicable, the presently reported expressions appear to be new.

The active medium is a collection of $`N`$ identical atoms, with the levels labeled $`|0`$, $`|1`$, $`|2`$, and $`|3`$ in increasing order of energy (see Fig. 1). Level separations are supposed to be large compared with $`k_\mathrm{B}T`$, where $`T`$ denotes the optical cavity temperature and $`k_\mathrm{B}`$ the Boltzmann constant, so that thermally-induced transitions are negligible. Note that the 3-level atoms schemes in Fig. 1(a) and 1(b) are *not* special cases of the 4-level scheme in Fig. 1(c). They are treated separately in appendices. Levels $`|1`$ and $`|2`$ are resonant with the field of a single-mode optical cavity. The active medium is supposed to be strongly homogeneously broadened so that atomic polarizations may be adiabatically eliminated, an approximation applicable also to YAG lasers, CO2 lasers and semiconductor lasers. Such lasers may operate in the steady-state regime, that is, pulsation and chaos do not normally occur.

The probability per unit time that an electronic transition from level $`|1`$ to level $`|2`$ occurs is taken as equal to $`m`$, and the probability of an electronic transition from $`|2`$ to $`|1`$ as $`m+1`$, where $`m`$ denotes the number of photons in the cavity (the qualification “per unit time” is henceforth omitted for the sake of brevity). This amounts to selecting a time unit whose typical value depends on the gain medium. Spontaneous decay from level $`|2`$ to level $`|1`$ is allowed with probability $`\gamma `$. This decay may be either non-radiative or involve radiation into other electromagnetic modes, besides the one of interest. Photons are absorbed with probability $`\alpha m`$, where $`\alpha `$ denotes a constant, the absorbing atoms residing most of the time in their ground state. These absorbing atoms model the transmission of light through mirrors with subsequent absorption by a detector. Provided detection is linear and reflectionless, it is immaterial whether absorption occurs inside or outside the optical cavity. For simplicity, internal absorption is neglected.

“Incoherent” pumping promotes electrons from level $`|0`$ to level $`|3`$ with probability $`P`$. When transitions from $`|0`$ to $`|3`$ and from $`|3`$ to $`|0`$ are both allowed with equal probabilities, the pumping process is called “coherent”, following an accepted terminology. We find it convenient to denote by $`\mathrm{}P`$ the $`|3|0`$ transition probability, with $`\mathrm{}=0`$ for incoherent pumping and $`\mathrm{}=1`$ for coherent pumping. Coherent pumping is physically realized by submitting the atoms to strong optical fields nearly resonant with the $`|0|3`$ transition. This pumping field may possibly originate from frequency-filtered thermal radiation and be highly multimode. Levels $`|0`$ and $`|3`$ need not be sharp. Instead, they may consist of narrow bands for improved coupling to broad-band pumps. One-way incoherent pumping would be appropriate to describe laser-diode pumps. In semiconductors, however, working levels spread into conduction and valence bands, so that the present model would not be adequate. Laser diodes have been treated previously on the basis of rate equations in Arnaud (2001). Spontaneous decay from level $`|3`$ to the upper working level $`|2`$ occurs with probability $`p_u`$, and spontaneous decay from the lower working level $`|1`$ to the ground level with probability $`p_d`$, where “u” stands for “up” and “d” for down. The relevant probabilities are schematized in Fig. 1.

One of the best-known laser-noise theory is probably that of Scully and Lamb (1967). Incoherent pumping is modeled by independent injection of 2-level atoms in the optical cavity. This model leads to a photo-count statistics which is, at best, Poissonian. More recently, Khazanov et al. (1990); Ralph and Savage (1991); Ritsch et al. (1991) considered the situation in which the pumping-levels populations may fluctuate. At first, it would seem that this may only increase the noise. It turns out, however, that population fluctuations are correlated in such a way that the output light fluctuations may be sub-Poissonian. It is difficult to pin point a simple intuitive explanation. It has been observed, however, that when lasers are pumped through a cascade of intermediate levels, pumping tends to be regular Ritsch and Zoller (1992); Briegel et al. (1996), a situation somewhat similar to laser-diodes high-resistance driving conditions.

Other means of generating sub-Poissonian light have been considered. Golubev and Sokolov (1984) were the first in 1984 to point out that lasers with non-fluctuating pumps should emit sub-Poissonian light. This conclusion has been verified experimentally by Machida et al. (1987) with the help of laser diodes driven by high-impedance electrical sources. The Scully and Lamb model has been generalized to account for regular atom injection Levien et al. (1993); Zhu (1993). Kolobov et al. (1993) made the interesting observation that the photodetection rate spectral density may be below the shot-noise level at *non-zero* Fourier frequencies in the case of *Poissonian* pumps. However, the photodetection rate spectral density remains at the shot-noise level at *zero* frequency. Accordingly, such lasers do not generate sub-Poissonian light in the sense defined earlier. It has been shown that 3-level lasers with coherent decay to the ground state Haake et al. (1993); Golubev (1996) generate sub-Poissonian light, and further that Raman lasers may generate sub-Poissonian light Ritsch et al. (1992). We will not consider here these more exotic configurations. A review is in Davidovitch (1996).

Rate equations treat the number of photons in the cavity as a classical random function of time. The light field is quantized as a result of matter quantization and conservation of energy, but not directly. Rate equations should be distinguished from semi-classical theories in which the optical field is driven by atomic dipole expectation values. The theory employed in this paper rests instead on the consideration of transition probabilities, as in the Loudon (1983) treatment of optical amplifier noise, for example. Every absorption event reacts on the number of light quanta in the optical cavity. Semi-classical theories are unable to explain sub-Poissonian light statistics because the light generation process and the light detection processes are considered separately. The expressions obtained from rate equations are found to coincide with Quantum Optics results when the number of atoms is large and transitions other than those relating to the atom-cavity interaction are incoherent. The laser is treated as a birth-death Markov process (see Sec. II). A Monte Carlo simulation gives the evolution of the number $`m`$ of photons in the cavity from which the Fano factor $`=\mathrm{var}(m)/m`$ is obtained. On the other hand, the instants $`t_k`$ when photons are being absorbed provide the spectral density of the photo-current, whose normalized value $`𝒮`$ is unity for Poisson processes. The normalized spectrum is denoted $`𝒮(\mathrm{\Omega })`$. The Fourier angular frequency $`\mathrm{\Omega }`$ is called for short: “frequency”.

When both the number of atoms and the pumping level increase, the computing time becomes prohibitively large because of the exponentially growing number of events to process. Analytical results are obtained by applying the weak-noise approximation. We first give in Sec. III.1 the steady-state photon number. The photo-current spectral density at zero frequency is evaluated in Sec. III.2, and the photo-current spectrum in Sec. III.3. The method of evaluation of the intra-cavity Fano factor is in Sec. III.4.

## II Monte Carlo method

The rate-equation model of a $`N`$–atoms singlemode laser straightforwardly leads to a master equation for the probability of having $`m`$ photons stored in the cavity at time $`t`$ Scully and Lamb (1967); Loudon (1983); Arnaud (2001). Alternatively, the laser evolution is modeled as a temporally homogeneous birth-death Markov process. In the steady-sate regime, rate of change and equilibrium probabilities of having $`m`$ and $`m+1`$ photons within the cavity are linked via a detailed balancing condition. This is a favorable condition for a Monte Carlo simulation Landau and Binder (2000) because every laser microstate belongs to a Markov chain and will thus occurs proportionally to its equilibrium probability when the number of step increases to $`\mathrm{}`$.

For the $`\mathrm{V}`$–type lasers considered in the present section (see Fig. 1(b)), rates of change $`W_j`$ are ascribed to the different kinds of events as given in Table 1. For example, the probability that an atom jumps from level $`|1`$ to level $`|3`$ during the elementary time interval $`[t,t+\delta t]`$ is $`W_2\delta t`$, where $`\delta t`$ is chosen small enough that this probability be much less than unity. Because atoms are coupled to one another only through the field, $`W_2`$ is proportional to the number $`n_1`$ of atoms in $`|1`$ at time $`t`$, and thus $`W_2=Pn_1`$, where the constant $`P`$ is proportional to the pump strength. If a jump does occur, $`n_1`$ is reduced by 1 while the number $`n_3`$ of atoms in level $`|3`$ is incremented by 1. If the initial value of $`n_3`$ is $`N`$ the event does not occur. Similar observations apply to the other jump probabilities. Notice that the coherent emission rate $`W_4`$ is proportional to $`m+1`$, following the Einstein prescription. This ensures that laser emission re-starts if extinction occurs. A key feature that distinguishes the present formulation from other rate-equation methods is that absorption of photons by the detector is included in the system description. Because detection is supposed to be linear, such events are taken to occur with a rate $`W_1=\alpha m`$, where $`\alpha `$ expresses detector absorption.

An efficient algorithm has actually been employed Gillespie (1992): Given that an event of any kind occurred at time $`\tau _k`$, the next-event time is

$$\tau _{k+1}=\tau _k+\frac{1}{_jW_j}\mathrm{ln}\left(\frac{1}{𝔯}\right),$$
(1)
where $`𝔯`$ is a random number uniformly distributed in the interval $`[0,1]`$. The probability that the event is of kind $`l`$ is equal to $`W_l/_jW_j`$. The Monte Carlo method was implemented on a desk computer. Only the total number of atoms in each state needs to be tracked (For fermions the numerical procedure is significantly more involved Chusseau and Arnaud (2001)). Within the whole time set $`\{\tau _k\}`$, we are mostly interested in the subset $`\{t_k\}`$ of photon-absorption events. It is straightforward in principle to evaluate the photo-detection noise. We also record $`m\left(t_k\right)`$ to evaluate the mean value and variance of $`m`$.

The intra-cavity Fano factor $``$ is represented in Fig. 2. Both coherent (Fig. 2(a)) and incoherent (Fig. 2(b)) pumping are considered. The analytical results in Sec. III.4, based on linearization agree well with the simulation. Notice that $``$ is below unity within some pumping range. There is good agreement with previous Quantum-Optics results Koganov and Shuker (2000).

The normalized spectral density $`𝒮(\mathrm{\Omega })`$ is represented in Fig. 3 for incoherent (a) and coherent (b) pumps, and two sets of parameter-values. For each Monte Carlo run, $`𝒮(\mathrm{\Omega })`$ is first evaluated from the $`\{t_k\}`$ list Chusseau and Arnaud (2001) and refined using a smoother power spectral density estimator Papoulis (1991); Press et al. (1992). Averaging over runs and concatenating neighboring frequencies produce the final data together with error bars at the 95% confidence level. There is fair agreement between Monte Carlo simulations and analytical formulas to be subsequently reported. Both predict sub-Poissonian photo-current statistics. Even with one billion photon-absorption events, Monte Carlo spectra exhibit large error bars. An analytical method is to be preferred when it exists. On the other hand, Monte Carlo simulations do not rely on linearization and provide a useful check.

## III Four-levels atoms lasers

### III.1 Steady-state

Let $`n_j`$, $`j=0,1,2,3`$, denote the number of atoms in state $`j`$, with

$$n_0+n_1+n_2+n_3=N.$$
(2)
The transition rate from $`|1|2`$ (stimulated absorption) is set equal to the number $`m`$ of photons in the cavity, a rule that defines a time unit. Pumping of an atom from level $`|0`$ to level $`|3`$ occurs with rate $`P`$. We allow for a transition rate $`\mathrm{}P`$ that an atom in level $`|3`$ decays back to level $`|0`$. Atoms in state $`|3`$ spontaneously decay to state $`|2`$ with rate $`p_u`$, while atoms in state $`|1`$ spontaneously decays to state $`|0`$ with rate $`p_d`$. These spontaneous decays are labelled “upper” and “lower” decay, respectively. Note that when $`1/p_u=0`$ the population of level $`|3`$ vanishes and any dependence on $`\mathrm{}`$ must therefore drop out. Photons are absorbed with rate $`\alpha m`$. Spontaneous decay from the upper to the lower working levels occurs with rate $`\gamma `$. Altogether, there are eight kinds of events that may occur in the course of time.

Let $`𝒥`$ denotes the net pumping rate, $``$ the net stimulated rate, $`𝒰`$ and $`𝒟`$ the upper and lower decay rates, $`𝒮`$ the spontaneous decay rate from the upper to the lower working levels, and $`𝒬`$ the photon absorption rate. The steady-state conditions then read

$`𝒥`$ $`=𝒰=𝒟=+𝒮,`$ (3a)
$`𝒬`$ $`=,`$ (3b)
where

$`𝒥`$ $`=Pn_0\mathrm{}Pn_3,`$ $`𝒰`$ $`=p_un_3,`$ (4a)
$``$ $`=(m+1)n_2mn_1,`$ $`𝒟`$ $`=p_dn_1,`$ (4b)
$`𝒮`$ $`=\gamma n_2,`$ $`𝒬`$ $`=\alpha m.`$ (4c)
Equations (2), (3) and (4) provide the steady-state atomic populations $`n_i`$ and photon number $`m`$. In particular

$$m=\frac{1}{2}\left(+\sqrt{^2+4𝒫𝒩}\right),$$
(5)
where

$`𝒩`$ $`={\displaystyle \frac{N}{\alpha }},`$ (6a)
$`𝒫`$ $`=\left[{\displaystyle \frac{1}{P}}+{\displaystyle \frac{2}{p_d}}+{\displaystyle \frac{1+\mathrm{}}{p_u}}\right]^1,`$ (6b)
$``$ $`=𝒫\left[𝒩1+{\displaystyle \frac{1}{p_d}}\left(1+\gamma 𝒩\gamma \right)\right]\gamma 1.`$ (6c)
For moderate pump powers above threshold, $`m`$ increases linearly with $`P`$. The intercept with the $`m=0`$ axis defines the threshold pump power. Mathematically, this amounts to first replacing in previous equations $`m+1`$ by $`m`$, and then setting $`m=0`$. The laser is found to oscillate for some pump rate provided

$$\gamma <\left(𝒩1\right)\left[\frac{𝒩+1}{p_d}+\frac{1+\mathrm{}}{p_u}\right]^1.$$
(7)
The limiting case is obtained at infinite pump power

$$\underset{P\mathrm{}}{lim}m=\frac{\left(𝒩1\right)\left(p_d\gamma \right)}{\left(1+\mathrm{}\right)\frac{p_d}{p_u}+2}\gamma .$$
(8)
Figure 4 shows how the steady-state photon number evolves as a function of the pumping parameter $`P`$ for various values of the spontaneous decay rate $`\gamma `$. Parameters are similar to Koganov and Shuker (2000) except $`p_u`$–values chosen here for a perfect superimposition of the coherent and incoherent plots. Moreover it corresponds to optimum pump-noise suppression as derived in Sec. III.2. At very high pump levels $`m`$ saturates like Eq. (8) because of ground-state population depletion.

### III.2 Zero-frequency noise

Our analytical results rest on a weak-noise approximation. Populations split into steady-state values and fluctuations. For example the instantaneous photon number $`m`$ is written as $`m+\mathrm{\Delta }m`$, where $`m`$ denotes the steady-state value. Rates split into steady-state values and fluctuations consisting of a deterministic function of the population fluctuations and uncorrelated Langevin “forces”. For example $`𝒥`$ split into $`J=𝒥`$ and $`\mathrm{\Delta }J`$. The latter is the sum of a deterministic function of the population fluctuations and a Langevin force $`j(t)`$ expressing the jump process randomness. Thus

$`𝒥`$ $`J+\mathrm{\Delta }J,`$ $`𝒰`$ $`U+\mathrm{\Delta }U,`$ (9a)
$``$ $`R+\mathrm{\Delta }R,`$ $`𝒟`$ $`D+\mathrm{\Delta }D,`$ (9b)
$`𝒮`$ $`S+\mathrm{\Delta }S,`$ $`𝒬`$ $`Q+\mathrm{\Delta }Q.`$ (9c)
where

$`\mathrm{\Delta }J=`$ $`P\mathrm{\Delta }n_0\mathrm{}P\mathrm{\Delta }n_3+j,`$ (10a)
$`\mathrm{\Delta }R=`$ $`\left(m+1\right)\mathrm{\Delta }n_2m\mathrm{\Delta }n_1`$
$`+\left(n_2n_1\right)\mathrm{\Delta }m+r,`$ (10b)
$`\mathrm{\Delta }S=`$ $`\gamma \mathrm{\Delta }n_2+s,`$ (10c)
$`\mathrm{\Delta }U=`$ $`p_u\mathrm{\Delta }n_3+u,`$ (10d)
$`\mathrm{\Delta }D=`$ $`p_d\mathrm{\Delta }n_1+d,`$ (10e)
$`\mathrm{\Delta }Q=`$ $`\alpha \mathrm{\Delta }m+q.`$ (10f)
A first-order variation of the expressions in Eq. (4) has been performed.

Conservation of the rates gives

$`\mathrm{\Delta }J`$ $`=\mathrm{\Delta }U=\mathrm{\Delta }D=\mathrm{\Delta }R+\mathrm{\Delta }S,`$ (11a)
$`\mathrm{\Delta }Q`$ $`=\mathrm{\Delta }R.`$ (11b)
Since the total number $`N`$ of atoms is constant we have

$$\mathrm{\Delta }n_0+\mathrm{\Delta }n_1+\mathrm{\Delta }n_2+\mathrm{\Delta }n_3=0.$$
(12)
Replacing atomic populations and photon numbers by their steady-state values in Sec. III.1, the above set of equations can be solved. In particular, $`\mathrm{\Delta }Q`$ is a linear combination of the Langevin forces

$$\mathrm{\Delta }Q=\underset{z\{j,d,u,q,r,s\}}{}c_zz,$$
(13)
where the $`c_z`$ are real coefficients that depend on the parameters $`N`$, $`P`$, $`\mathrm{}`$, $`p_u`$, $`p_d`$, $`\gamma `$ and $`\alpha `$. The detailed expressions, too lengthy to be given here in their general form, are conveniently handled using symbolic calculations.

The normalized zero-frequency photo-current spectral density is of the form

$$𝒮=\frac{1}{\alpha m}\underset{z\{j,d,u,r,s\}}{}c_z^2\sigma _z,$$
(14)
where $`\sigma _z`$ denotes the spectral density value of the Langevin noise source $`z`$, equal to average rates.

$`\sigma _j`$ $`=Pn_0+\mathrm{}Pn_3,`$ $`\sigma _u`$ $`=p_un_3,`$ (15a)
$`\sigma _r`$ $`=\left(m+1\right)n_2+mn_1,`$ $`\sigma _d`$ $`=p_dn_1,`$ (15b)
$`\sigma _s`$ $`=\gamma n_2,`$ $`\sigma _q`$ $`=\alpha m.`$ (15c)
When these expressions are introduced in Eq. (14) an analytical expression of $`𝒮`$ is obtained. Three special cases are considered below, (a) $`\gamma =0`$ and $`m`$ large compared with unity, (b) $`N\alpha `$, (c) $`\gamma =0`$ and $`N\alpha `$.

(a) If spontaneous decay is negligible ($`\gamma =0`$), Eq. (14) yields

$$\begin{array}{c}𝒮=1+\frac{2}{\left(𝒩1\right)^2}+\frac{8𝒫^2}{p_d^2}+\frac{\left(64𝒩\right)𝒫}{\left(𝒩1\right)p_d}\hfill \\ \hfill +\frac{2\left(1+\mathrm{}\right)𝒫^2}{p_u^2}+\frac{2𝒫\left(2𝒫p_d\right)}{p_dp_u},\end{array}$$
(16)
where $`𝒫`$ and $`𝒩`$ are defined in Eq. (6). The normalized spectral density is unity at low and high pumping levels. For some constant $`𝒩`$ value, $`𝒮`$ reaches its minimum valued

$$𝒮_{min}=\frac{2𝒩\left(𝒩1\right)+11+\mathrm{}\left(4𝒩\left(𝒩1\right)+15\right)}{2\left(3+4\mathrm{}\right)\left(𝒩1\right)^2},$$
(17)
when

$`{\displaystyle \frac{P}{p_u}}`$ $`={\displaystyle \frac{1}{1+\mathrm{}}},`$ $`{\displaystyle \frac{P}{p_d}}`$ $`={\displaystyle \frac{𝒩\left(1+2\mathrm{}\right)\left(2+3\mathrm{}\right)}{\left(1+\mathrm{}\right)\left(2𝒩1\right)}}.`$ (18)
(b) When $`N\alpha `$, Eq. (14) yields

$$\begin{array}{c}𝒮=1+\frac{2\gamma }{p_d\gamma }\frac{4𝒫+2\gamma }{p_d}+\frac{8𝒫\left(𝒫+\gamma \right)}{p_d^2}\frac{8𝒫^2\gamma }{p_d^3}\hfill \\ \hfill +\frac{2𝒫\left(\gamma p_d\right)\left(p_d2𝒫\right)}{p_d^2p_u}\frac{2\left(1+\mathrm{}\right)𝒫^2\left(\gamma p_d\right)}{p_dp_u^2}.\end{array}$$
(19)
Figure 5 gives the normalized zero-frequency photo-current spectral density $`𝒮(P/p_d,\gamma /p_d)`$ in the form of contour plots, selecting $`P/p_u=2P/p_d`$ for the case of incoherent pumping and $`P/p_u=\frac{2}{3}P/p_d`$ for the case of coherent pumping. The darker the area, the lower is the spectral density. Since dark areas are wider in Fig. 5(a) than in Fig. 5(b), incoherent pumping is to be preferred. Figure 5 shows that sub-Poissonian light generation by optically pumped 4-level atom lasers is robust against spontaneous decay and pumping conditions. The optimum conditions (darkest areas) are defined in Eq. (18). But small departures from these conditions do not increase much the noise.

(c) If both $`\gamma =0`$ and $`N\alpha `$, the spectral density obtained either by setting $`𝒩=\mathrm{}`$ in Eq. (16) or $`\gamma =0`$ in Eq. (19), reads

$$𝒮=1\frac{2Pp_dp_u\left(p_d+2\left(P+P\mathrm{}+p_u\right)\right)}{\left(2Pp_u+p_d\left(P+P\mathrm{}+p_u\right)\right)^2},$$
(20)
an expression that coincides with Eq. (4) of Ritsch et al. (1991). The absolute minimum value and corresponding conditions are obtained from Eqs. (17) and (18)

$`𝒮_{min}`$ $`={\displaystyle \frac{1+2\mathrm{}}{3+4\mathrm{}}},`$ $`{\displaystyle \frac{P}{p_u}}`$ $`={\displaystyle \frac{1}{1+\mathrm{}}},`$ $`{\displaystyle \frac{P}{p_d}}`$ $`={\displaystyle \frac{1+2\mathrm{}}{2\left(1+\mathrm{}\right)}}.`$ (21)
For incoherent pumping, $`\mathrm{}=0`$, we have therefore $`𝒮_{min}=1/3`$ when $`p_u=P`$ and $`p_d=2P`$. For coherent pumping, $`\mathrm{}=1`$, we have $`𝒮_{min}=3/7`$ when $`p_u=2P`$ and $`p_d=\frac{4}{3}P`$.

Table 2 gives the minimum spectral density values achievable with optically pumped three- and 4-level lasers as obtained from Eqs. (21), (38) and (32). Under the very special condition of negligible spontaneous decay and $`N\alpha `$, the intra-cavity Fano factor depends linearly of the zero-frequency normalized photo-current spectral density Golubev and Sokolov (1984), $`𝒮=21`$. Such relation does not hold in the general situation where Fano factors follow from the formulation in Sec. III.4.

$`\mathrm{V}`$–type incoherently-pumped lasers were treated earlier by Khazanov et al. (1990). Ralph and Savage (1991, 1993) extended the analysis to incoherently pumped $`\mathrm{\Lambda }`$–type lasers and 4-level atom lasers. Ritsch et al. (1991) gave a description of 4-level lasers for the two pumping schemes. These previous results are exactly recovered from the present rate-equations method.

Formulas derived in the present paper for coherently-pumped 3-level lasers, however, appear to be new. Because the present method involves simple algebra the conditions under which the photo-current spectral density is minimum are easily obtained. To our knowledge, all the results presented in this paper when $`\gamma `$ does not vanish, or $`N`$ is not much larger than $`\alpha `$, are new.

### III.3 Photo-detection spectrum

The steady-state relations were given in Eq. (11). At some Fourier frequency $`\mathrm{\Omega }`$ the generalized rate equations read Arnaud (1995)

$`i\mathrm{\Omega }\mathrm{\Delta }m`$ $`=\mathrm{\Delta }R\mathrm{\Delta }Q,`$ (22a)
$`i\mathrm{\Omega }\mathrm{\Delta }n_0`$ $`=\mathrm{\Delta }D\mathrm{\Delta }J,`$ (22b)
$`i\mathrm{\Omega }\mathrm{\Delta }n_1`$ $`=\mathrm{\Delta }R+\mathrm{\Delta }S\mathrm{\Delta }D,`$ (22c)
$`i\mathrm{\Omega }\mathrm{\Delta }n_2`$ $`=\mathrm{\Delta }U\mathrm{\Delta }R\mathrm{\Delta }S,`$ (22d)
$`i\mathrm{\Omega }\mathrm{\Delta }n_3`$ $`=\mathrm{\Delta }J\mathrm{\Delta }U,`$ (22e)
Similar to Sec. III.2, equations Eqs. (9), (10), (12) and (22) are solved for $`\mathrm{\Delta }Q`$. The formula for the light spectral density $`𝒮`$ is the same as Eq. (14) except that the coefficients $`\stackrel{~}{c}_z`$ are complex and frequency dependent

$$𝒮(\mathrm{\Omega })=\frac{1}{\alpha m}\underset{z\{j,d,u,r,s\}}{}\stackrel{~}{c}_z(\mathrm{\Omega })\stackrel{~}{c}_z^{}(\mathrm{\Omega })\sigma _z,$$
(23)
where the Langevin “forces” $`\sigma _z`$ are still given in Eq. (15).

After rearranging, Eq. (23) gives the spectral density in the form

$$𝒮(\mathrm{\Omega })=1+\frac{a_3\mathrm{\Omega }^6+a_2\mathrm{\Omega }^4+a_1\mathrm{\Omega }^2+a_0}{\mathrm{\Omega }^8+b_3\mathrm{\Omega }^6+b_2\mathrm{\Omega }^4+b_1\mathrm{\Omega }^2+b_0},$$
(24)
where the coefficients $`a_i`$ and $`b_i`$ are real. The form in Eq. (24) ensures that $`𝒮(\mathrm{\Omega })`$ tends to unity (shot-noise level) at high frequencies.

Figure 6 shows that $`𝒮`$ reaches its minimum value at $`\mathrm{\Omega }=0`$. When spontaneous decay from the upper working level may be neglected, light is always sub-Poissonian and the lowest $`𝒮`$–value occurs when $`P/p_d=1/2`$. Spontaneous decay from the upper working level is inconsequential until $`\gamma /p_d\mathrm{3\hspace{0.17em}10}^2`$. The light statistics ceases to be sub-Poissonian when $`P/p_d>0.3`$.

### III.4 Fano factor

The intra-cavity photon statistics is characterized by the Fano factor $`=\mathrm{\Delta }m^2/m`$, where $`m`$ denotes as before the number of photons in the cavity. The variance of $`m`$ equals the integration over frequency of $`𝒮_{\mathrm{\Delta }m}(\mathrm{\Omega })`$, the spectral density of $`\mathrm{\Delta }m`$. Similarly to Sec. III.3, $`𝒮_{\mathrm{\Delta }m}(\mathrm{\Omega })`$ obtains by solving for $`\mathrm{\Delta }m`$ instead of $`\mathrm{\Delta }Q`$ and its general form is

$$𝒮_{\mathrm{\Delta }m}(\mathrm{\Omega })=\frac{a_3^{}\mathrm{\Omega }^6+a_2^{}\mathrm{\Omega }^4+a_1^{}\mathrm{\Omega }^2+a_0^{}}{\mathrm{\Omega }^8+b_3^{}\mathrm{\Omega }^6+b_2^{}\mathrm{\Omega }^4+b_1^{}\mathrm{\Omega }^2+b_0^{}},$$
(25)
where $`a_i^{}`$ and $`b_i^{}`$ are real coefficients. The Fano factor is thus

$$=_{\mathrm{}}^{\mathrm{}}𝒮_{\mathrm{\Delta }m}(\mathrm{\Omega })\frac{\mathrm{d}\mathrm{\Omega }}{2\pi }.$$
(26)
The Fano factor for incoherently pumped 4-level atoms lasers is represented as a function of pump and spontaneous decay rates in Fig. 7.

At small pumping rates, and when the spontaneous decay rate $`\gamma `$ is large, the photon statistic is essentially that of amplified thermal light and the Fano factor $`=m+1`$. The inset in Fig. 7 shows that this simple result holds indeed for $`\gamma =632`$, excepts for the highest values of $`P`$, where a slight reduction of $``$ occurs.

When $`\gamma <p_d`$, the Fano factor exhibits a peak just below threshold Jones et al. (1999), and decreases to $`1`$ until the ground level gets significantly depleted. Self-quenching Mu and Savage (1992); Koganov and Shuker (2000) is not observed here because atomic polarizations have been adiabatically eliminated. Otherwise, our results fully agree with those reported by Koganov and Shuker (2000).

## IV Conclusion

We have considered optically-pumped 4-level and 3-level atom lasers in resonant single-mode cavities. The light statistics has been obtained from a simple rate-equation approach, using both a Monte Carlo simulation and an analytical method based on linearization. The emitted light may be sub-Poissonian, as was previously observed by many authors. Whenever a comparison can be made, exact agreement is reached with the previous Quantum Optics results. In the case of coherently-pumped 3-level atoms lasers our results appear to be new. Note that 3-level atoms lasers are not special cases of the 4-level scheme. When the assumptions of negligible spontaneous decay and large atom numbers are not made, the results presented in this paper for the internal and external field statistics appear to be new.

For practical reasons, Monte Carlo simulations were restricted to $`N1000`$ atoms. Because the analytical formulas, obtained through the use of symbolic calculus, are lengthy they were not written down in the paper. However, they were employed to determine the conditions under which the spectral density of the photo-current reaches its minimum value. For example, we found that when spontaneous decay from the upper working level may be neglected, 3-level atoms lasers may deliver light with fluctuations at half the shot-noise level. 4-level atoms lasers may deliver light with fluctuations at one third of the shot-noise level. The photo-current noise decreases further and tends to zero, under ideal conditions, when the number of levels becomes large Briegel et al. (1996).

## Appendix A  $`\mathrm{\Lambda }`$–type lasers

The number $`m`$ of photons in the steady-state is still

$$m=\frac{1}{2}\left(+\sqrt{^2+4𝒫𝒩}\right),$$

with

$`𝒩`$ $`={\displaystyle \frac{N}{\alpha }},`$ (27a)
$`𝒫`$ $`=\left[{\displaystyle \frac{1}{P}}+{\displaystyle \frac{2+\mathrm{}}{p_d}}\right]^1,`$ (27b)
$``$ $`=𝒫\left[𝒩1\mathrm{}+{\displaystyle \frac{1}{p_d}}\left(\left(1+\mathrm{}\right)\left(1+\gamma \right)𝒩\gamma \right)\right]`$
$`1\gamma .`$ (27c)
Lasing may occur if

$$\gamma <\frac{𝒩1\mathrm{}}{𝒩+1}p_d.$$
(28)
(a) When spontaneous decay may be neglected, $`\gamma =0`$, the zero-frequency normalized photo-current spectral density reads

$$\begin{array}{c}𝒮=1+\frac{2\left(1+\mathrm{}+\mathrm{}𝒩\right)}{\left(1+\mathrm{}𝒩\right)^2}+\frac{4\left(2+\mathrm{}\right)𝒫^2}{p_d^2}\hfill \\ \hfill \frac{2\left(3+2\mathrm{}2𝒩\right)𝒫}{\left(1+\mathrm{}𝒩\right)p_d}.\end{array}$$
(29)
Note that when $`\mathrm{}=0`$ this expression may be obtained by setting $`1/p_u=0`$ in the 4-level expression. This is not, however, a valid procedure in general.

(b) When $`N\alpha `$

$$\begin{array}{c}𝒮=1+\frac{2\gamma }{p_d\gamma }\frac{2\left(2𝒫+\gamma \right)}{p_d}\hfill \\ \hfill +\frac{2𝒫\left(2\left(2+\mathrm{}\right)𝒫+\left(4+\mathrm{}\right)\gamma \right)}{p_d^2}\frac{4\left(2+\mathrm{}\right)𝒫^2\gamma }{p_d^3}.\end{array}$$
(30)
(c) When $`\gamma =0`$ and $`N\alpha `$

$$𝒮=1\frac{4Pp_d}{\left(P\left(2+\mathrm{}\right)+p_d\right)^2}$$
(31)
The minimum value of $`𝒮`$ and corresponding value of $`P`$ are

$`𝒮_{min}`$ $`={\displaystyle \frac{1+\mathrm{}}{2+\mathrm{}}},`$ $`P`$ $`={\displaystyle \frac{p_d}{2+\mathrm{}}}.`$ (32)
## Appendix B $`\mathrm{V}`$–type lasers

The number $`m`$ of photons in the steady-state is still

$$m=\frac{1}{2}\left(+\sqrt{^2+4𝒫𝒩}\right),$$

with

$`𝒩`$ $`={\displaystyle \frac{N}{2\alpha }},`$ (33a)
$`𝒫`$ $`=\left[{\displaystyle \frac{1}{P}}+{\displaystyle \frac{1+2\mathrm{}}{2p_u}}\right]^1,`$ (33b)
$``$ $`={\displaystyle \frac{𝒫}{4p_u}}\left[\left(2𝒩1\right)\left(\gamma +2p_u\right)1\right]`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\gamma }{2}}𝒩\gamma .`$ (33c)
Lasing may occur if

$$\gamma <\frac{2𝒩1}{1+\mathrm{}+2\mathrm{}𝒩}p_u.$$
(34)
(a) When spontaneous decay may be neglected, $`\gamma =0`$, the zero-frequency normalized photo-current spectral density reads

$$\begin{array}{c}𝒮=1+\frac{2𝒩+1}{\left(2𝒩1\right)^2}+\frac{\left(1+2\mathrm{}\right)𝒫^2}{2p_u^2}\hfill \\ \hfill \frac{\left(4𝒩1\right)𝒫}{2\left(2𝒩1\right)p_u}.\end{array}$$
(35)
(b) When $`N\alpha `$

$$\begin{array}{c}𝒮=1+\frac{2\gamma }{𝒫\gamma }\frac{𝒫}{p_u}+\frac{𝒫\left(𝒫\left(1+2\mathrm{}\right)\gamma \right)}{2p_u^2}\hfill \\ \hfill +\frac{\left(1+2\mathrm{}\right)𝒫^2\gamma }{4p_u^3}\frac{2𝒫^2\gamma }{\left(𝒫\gamma \right)\left(𝒫\gamma +2\left(𝒫\gamma \right)p_u\right)}.\end{array}$$
(36)
(c) When $`\gamma =0`$ and $`N\alpha `$

$$𝒮=1\frac{4Pp_u}{\left(P\left(1+2\mathrm{}\right)+2p_u\right)^2}$$
(37)
The minimum value of $`𝒮`$ and corresponding value of $`P`$ are

$`𝒮_{min}`$ $`={\displaystyle \frac{1+4\mathrm{}}{2+4\mathrm{}}},`$ $`P`$ $`={\displaystyle \frac{2p_u}{1+2\mathrm{}}}.`$ (38)
###### Acknowledgements.

This work was supported by the STISS Department of Université Montpellier II and by CNRS under the JemSTIC Program.