text:
# Two Simple Approaches to Sol-Gel Transition

## 1 INTRODUCTION

Within the phenomenological framework the basic distinction between a liquid and a glass consists in character of the relaxation law of shear components of the elastic stress: if in an ideal glass they are kept infinitely long, in a liquid such relaxation proceeds in the finite time $`\tau =\eta /G`$, where $`\eta `$ is the dynamical shear viscosity and $`G`$ is the shear modulus . In naive manner it is possible to assume, that a glass transition is caused purely by kinetic effect of a liquid freezing, for which viscosity $`\eta `$ gets infinite value for a finite value of the shear modulus $`G`$ . However in the course of usual second-order phase transition, where infinite increase of the $`\tau `$ at critical point is also observed the situation is reverse. Really, proceeding from the viscoelastic liquid to a general case, one has $`\tau =\chi /\gamma `$, where $`\chi `$ is the generalized susceptibility and $`\gamma `$ is the kinetic coefficient (in the case under consideration they are reduced to the quantities $`G^1,\eta ^1`$, respectively) . The infinite increase of a susceptibility $`\chi `$ occurs and a kinetic coefficient $`\gamma `$ does not manifest any peculiarity at the phase transition. In our case, this is equivalent to the situation when the shear modulus $`G`$ is approaching zero value at a fixed viscosity $`\eta `$. This situation corresponds to a viscoelastic transition .

Usually under a glass transition, the transition into a glassy state which occurs during fast cooling of a liquid is understood. However, an analogous type of transition occurs during polymer network formation (polymer gelation) from monomers or linear macromolecules (sol) by means of a chemical reaction at constant temperature. This transition is referred to as the sol-gel transition. In such a case, thermodynamic peculiarities are observed as, e.g., appearance of the shear modulus when the system reaches a critical point.111 Of course, when the temperature of the polymer network formed is decreased the kinetic glass transition in above sense occurs, as well. The first theoretical representations of such type transitions have been elaborated in classical works by Flory . Modern description of the glassy state of the macromolecule networks (see ) is based on pioneering contributions to the theory of soft condensed matter: the Deam-Edwards theory  of a cross-linked tangled macromolecule, and the Edwards-Anderson theory  of spin glass.

The aim of this paper is to elaborate the simplest theoretical scheme so that the transition occurring during polymer gelation can be explained on the base of obvious combinatoric arguments (Section 2) as well as within framework of a synergetic theory (Section 3). In the former case, the critical point and behaviour of the system behind the critical point are obtained using information about the cross-linking properties of the constituent parts of the system (monomers) and extent of the cross-linking process, only. In Section 3, the starting point is that a sol-gel transition is ensured by self-organization of the elastic fields of stress and strain shear components and the number of cross-links. Both the combinatoric approach and the synergetic one allow to obtain anomalous dependence of the shear modulus on the number of cross-links.

## 2 Combinatoric approach

The gelation process has been described first by Flory and Stockmayer. They used a laborious approach based on combinatoric considerations of the most probable composition of the system. A much more effective variant of this approach, exploiting theory of branching processes, has been developed later by Dobson and Gordon . This approach allows also an analysis of structural details of the system behind the critical point, see, e.g. .

As the simplest case let us consider the gelation process in the system which at the beginning consists of a large number of monomers $`𝒩`$, wearing $`f`$ functional terminal groups of independent reactivity. Monomers react mutually and irreversibly via their terminal groups and if $`f>2`$, branched molecules of increasing size and complexity are formed progressively in the system. Extent of the reaction is described by the conversion of terminal groups $`\alpha `$, which is defined as the ratio of the number of groups consumed by the reaction at given time and the starting number of groups. Eventually, when a critical conversion $`\alpha _\mathrm{c}`$ is achieved a molecule of macroscopic dimensions (gel) appears in the system , . Gel contains cycles of long sequences of linked monomers and, consequently, attains elastic properties. Flory  has shown that shear modulus of gel $`G`$ is proportional to its “cycle rank” $`\xi `$ defined as the number of “superfluous” links formed in gel which can be cut without breaking integrity of gel. In other words, cycle rank is the number of cuts needed for elimination of all cycles from gel.

The key role in the approach is played by the extinction probability $`v`$, which is a probability that a link formed in the process has just a finite continuation. It can be shown easily  that for the system considered here the extinction probability can be obtained as a root of the equation

$$v=(1\alpha +\alpha v)^{f1}$$
(1)
satisfying the condition $`0v1`$. Eq.(1) expresses fact that a given link has a finite continuation in a direction only if $`f1`$ functional groups remaining on the monomer connected by the link are either unreacted (with probability $`1\alpha `$) or reacted giving rise to links with finite continuation only (with probability $`\alpha v`$). Below the critical conversion, the only solution of Eq.(1) is $`v=1`$, i.e., only molecules of finite size are formed in the system. However, behind the critical conversion monomers can be found either in sol or in gel: $`v<1`$. As a measure of the “distance” from the critical conversion let us introduce the parameter $`ϵ`$ defined as

$$ϵ\alpha \alpha _\mathrm{c}$$
(2)
and expand the extinction probability in a series:

$$v=1+A_1ϵ+A_2ϵ^2+\mathrm{}$$
(3)
Substituting Eqs.(2) and (3) into Eq.(1), expressions for the critical conversion and the parameters $`A_i`$ are obtained as

$`\alpha _\mathrm{c}`$ $`=`$ $`{\displaystyle \frac{1}{f1}},`$ (4)
$`A_1`$ $`=`$ $`2{\displaystyle \frac{(f1)^2}{f2}},`$ (5)
$`A_2`$ $`=`$ $`{\displaystyle \frac{4}{3}}{\displaystyle \frac{(2f3)(f1)^3}{(f2)^2}}.`$ (6)
The expression for the cycle rank can be found in the following way. Obviously, to join $`𝒩`$ monomers into a cycle-free structure, $`𝒩1`$ links are necessary. By definition, the cycle rank $`\xi `$ is the number of “superfluous” links formed in gel, i.e., the difference between the number $`\mathrm{N}_\mathrm{G}`$ of all links formed in gel and the number $`𝒩_\mathrm{G}1`$ of links sufficient to join together $`𝒩_\mathrm{G}`$ monomers of gel:

$$\xi \mathrm{N}_\mathrm{G}(𝒩_\mathrm{G}1)\mathrm{N}_\mathrm{G}𝒩_\mathrm{G}$$
(7)
as $`\mathrm{N}_\mathrm{G}`$, $`𝒩1`$. Correspondingly, the number of links formed in gel is the difference between the numbers of links formed in the total system and in sol, i.e.,

$$\mathrm{N}_\mathrm{G}=\frac{1}{2}𝒩\alpha f(1v^2)$$
(8)
as any link in sol has to have finite continuations in two directions. On the other hand, the number $`𝒩_\mathrm{G}`$ of monomers incorporated in gel is the difference between the total number of monomers and the number of monomers in sol which is made of monomers with links of only finite continuation:

$$𝒩_\mathrm{G}=𝒩𝒩(1\alpha +\alpha v)^f.$$
(9)
By virtue of Eqs.(7) – (9), one gets for cycle rank of the system considered

$$\frac{\xi }{𝒩}=(1\alpha +\alpha v)^f+\frac{1}{2}\alpha f(1v^2)1.$$
(10)
Finally, substituting Eqs.(3),(5) and (6) into the formula (10), one gets the necessary expansion:

$$\frac{\xi }{𝒩}=\frac{2}{3}\frac{(f1)^4}{(f2)^2}fϵ^3+O(ϵ^4).$$
(11)
Respectively, the weight fraction of gel $`w_\mathrm{G}𝒩_\mathrm{G}/𝒩`$ is determined by Eq.(9) to read

$$w_\mathrm{G}=2\frac{(f1)^2}{f2}ϵ+O(ϵ^2).$$
(12)
So, if the critical exponent of the gel weight fraction is equal 1 as usual, this for the cycle rank is anomalous large being equal 3. It is worthwhile to note that combinatoric approach is based exclusively on the information about the functionality of the monomers and extent of the chemical reaction between the monomers in the system considered.

## 3 Synergetic approach

Now, let us consider the polymer network as a viscoelastic continuum matter that is characterized by the shear modulus $`G`$ and the shear viscosity $`\eta `$. Process of polymer gelation is determined by the number of the cross-links $`N`$, which value is different from a stationary magnitude $`N_0`$ at a time $`t`$. Therein, an elastic state of the polymer is defined by the shear component of the proper (internal) values of deformation $`\epsilon (t)`$ and stress $`\sigma (t)`$. The keypoint is that these values are not reduced to the external elastic deformation $`e1`$ and stress $`\sigma _eG`$, in particular they can get large magnitudes $`\epsilon 1`$, $`\sigma G`$.

Our consideration of evolution of the elastic continuum with the internal structure is stated on the phenomenological equations by Maxwell-Kelvin

$$\frac{d\epsilon }{dt}=\frac{\epsilon }{\tau }+\frac{\sigma }{\eta },$$
(13)
$$\frac{d\sigma }{dt}=\frac{\sigma }{\tau _\sigma }+g_\sigma \epsilon N.$$
(14)
Here we introduce a macroscopic relaxation time $`\tau `$ for the strain and a microscopic one $`\tau _\sigma `$ for the stress, as well as a constant $`g_\sigma >0`$ of the positive feedback between the deformation $`\epsilon `$ and the number of cross-links $`N`$. Within the microscopic interval $`t\tau _\sigma `$, steady-state condition $`d\sigma /dt=0`$ in Eq.(14) leads to the Hooke law with the microscopic shear modulus

$$G_\sigma \tau _\sigma g_\sigma N$$
(15)
being determined by the number of cross-links $`N`$. Respectively, within a macroscopic interval $`t\tau `$ Eq.(13) gives the magnitude $`G\eta /\tau `$ that is characteristic for the usual modulus of the viscoelastic matter. Lastly, a variation rate $`dN/dt`$ of the internal degree of freedom is supposed to be determined by the equation

$$\frac{dN}{dt}=\frac{N_0N}{\tau _N}g_N\sigma \epsilon $$
(16)
where $`\tau _N`$ is a mesoscopic relaxation time, $`g_N>0`$ is constant of negative feedback between the deformation $`\epsilon `$ and the stress $`\sigma `$. Within a mesoscopic interval $`\tau _Nt\tau `$, Eq.(16) determines a steady-state value

$$N=N_0\tau _Ng_N\sigma \epsilon $$
(17)
that is smaller than the magnitude $`N_0`$ fixed by external conditions due to the fact that the elastic energy is proportional to the product $`\sigma \epsilon `$.

System of Eqs.(13), (14) and (16) is known in synergetics  as the Lorenz system where the deformation $`\epsilon `$, the stress $`\sigma `$ and the number of cross-links $`N`$ play roles of an order parameter, a conjugate field and a control parameter, respectively. It is very important for following considerations that the relation between micro-, meso- and macroscopic values of the relaxation times

$$\tau _\sigma ,\tau _N\tau $$
(18)
is satisfied. Due to this condition the evolution of the quantities $`\sigma `$, $`N`$ turns out to be subordinated to the long-time variation of $`\epsilon `$. A peculiarity of the Lorenz system consists in linear character of the equation (13) for the order parameter $`\epsilon `$ and in non-linearity of equations (14), (16) for the conjugate field $`\sigma `$ and the control parameter $`N`$. The negative nature of non-linearity in Eq.(16) means a decrease of the number $`N`$ of cross-links. Evidently, this fact reflects Le Chatelier principle. A non-linear term in Eq.(14) for a field $`\sigma `$ describes the positive feedback causing the system self-organization.

Expressions (13), (14) and (16) form the complete system of equations determining the polymer cross-linking behaviour. Because of a slow evolution, the order parameter $`\epsilon (t)`$ subordinates variations of quantities $`\sigma (t)`$, $`N(t)`$, so that one can take $`d\sigma /dt=dN/dt=0`$ within the framework of the adiabatic approximation . Then $`N`$, $`\sigma `$ are expressed in terms of $`\epsilon `$ by the equations:

$$N=\frac{N_0}{1+\epsilon ^2/\epsilon _m^2},\epsilon _m^2\tau _\sigma \tau _Ng_\sigma g_N;$$
(19)
$$\sigma =G_0\frac{\epsilon }{1+\epsilon ^2/\epsilon _m^2},G_0\tau _\sigma g_\sigma N_0.$$
(20)
In accordance with Eq.(19) the number of cross-links $`N`$ decreases monotoneously with increase of the strain $`\epsilon `$ from the value $`N_0`$ at $`\epsilon =0`$ to $`N_0/2`$ at $`\epsilon =\epsilon _m`$.222Obviously this decrease is caused by the negative feedback in Eq.(16), that is reflection of Le Chatelier principle for analyzes problem. Actually, a liquid self-organization, resulting in a sol-gel transition, is caused by the positive feedback between the strain $`\epsilon `$ and the number of cross-links $`N`$ in Eq.(14). Hence, the increase of the number of cross-links should intensify the self-organization effect. However, according to Eq.(19) system behaves in such way that the consequence of self-organization, i.e., growth of the elastic strain, leads to decrease of its origin – the number of cross-links. In Eq.(20) the elastic stress in terms of the strain has the linear form of the Hooke law at $`\epsilon \epsilon _m`$ with the effective shear modulus $`G_0`$. Then, at $`\epsilon =\epsilon _m`$ the function $`\sigma (\epsilon )`$ has a maximum and at $`\epsilon >\epsilon _m`$ it decreases which is physically meaningless. Thus, the constant $`\epsilon _m`$, defined by the second equation (19), has the meaning of the maximum achievable strain.

Substituting Eq.(20) in Eq.(13) we find the equation describing evolution of a system in the course of the sol-gel transition:

$$\frac{d\epsilon }{dt}=\gamma \frac{E}{\epsilon },\gamma \frac{\epsilon _m^2}{\tau TN_0}$$
(21)
where constant $`\gamma `$ plays a role of the kinetic coefficient. Behaviour of the system under consideration is determined by the dependence $`E(\epsilon )`$ of the elastic energy on the strain:

$$E\frac{TN_0}{2}\left[\frac{\epsilon ^2}{\epsilon _m^2}\frac{N_0}{N_c}\mathrm{ln}\left(1+\frac{\epsilon ^2}{\epsilon _m^2}\right)\right]$$
(22)
where the characteristic value of the number of cross-links is introduced

$$N_c\frac{\eta }{\tau _\sigma g_\sigma }.$$
(23)
At $`N_0N_c`$ dependence (22) is monotoneously increasing with a minimum at the point $`\epsilon =0`$. It means that in the stationary state ($`\dot{\epsilon }=0`$) the elastic strain is not spontaneous. Thus, a liquid state is realized, in which the strain caused by the external stress relaxes during the time

$$\tau _{ef}=\tau \left(1N_0/N_c\right)^1.$$
(24)
The relaxation time increases infinitely when the number of cross-links $`N_0`$ reaches the critical value $`N_c`$ and at $`N_0>N_c`$ the system undergoes a sol-gel transition. In the gel state the multiplier 1/2 appears in the dependence (24), and the minimum of the elastic energy corresponds to the elastic strain

$$\frac{\epsilon ^2}{\epsilon _m^2}=\frac{N_c}{N_0}\frac{N_0N_c}{N_c}\frac{ϵ}{1+ϵ}$$
(25)
where we introduce the distance from the critical value $`N_c`$

$$ϵ\frac{N_0N_c}{N_c}$$
(26)
being equivalent to the definition (2). Inserting Eq.(24) into the dependence (22), we obtain the elastic energy of the steady-state:

$$E_0E(\epsilon _0)=\frac{TN_0}{4}ϵ^2+O(ϵ^3).$$
(27)
As would be expected, this energy is proportional to the second power of the parameter (26) and is negative in nature (the latter means the energy benefit of the gel state in comparison with the liquid state).

Taking into account that the glassy state is determined by density of the localized monomers, let us find now the shear modulus of the appeared gel state. It is principally important in our considerations that the gel state is determined by the value of the Deam-Edwards parameter of localization $`\omega `$  which is supposed to be proportional to the square of the proper strain $`\epsilon _0`$ of the matter. Under the condition of the appearance of the elastic strains $`e`$, a generalized Deam-Edwards parameter $`\omega (e)`$ has to be considered which is related to the condition $`\omega (e=0)\omega `$. Then, expanding the function $`\omega (e)`$ into a series and keeping the first two terms only, one obtains

$$\omega (e)\omega (1+e^2),\omega \epsilon _0^2,e1.$$
(28)
By virtue of the parity condition $`\omega (e)=\omega (e)`$, this expansion does not contain a linear term.333 Otherwise, the usual squared dependence (31) will not be obtained. Because of that the total strain $`\epsilon \sqrt{\omega (e)}`$, internal one $`\epsilon _0\sqrt{\omega }`$ and elastic one $`e`$ are connected by the following relation:

$$\epsilon ^2=\epsilon _0^2(1+e^2),e1.$$
(29)
The key point is that Eq.(29) supposes the additivity rule holds not for quantities $`\epsilon `$, $`\epsilon _0`$, $`e`$ themselves, but for their squares. The physical reason for such a situation is that the system under consideration is random in character and described by a symmetrical distribution function. Therefore, all odd-power moments vanish identically and making use of the additivity rule (29) for variances follows.

It is worthwhile to note the seeming contradiction between above relations for the localization parameter $`\omega `$ and the Deam-Edwards results . Obvious reason consists in that the formers are obtained within framework of the mean-field theory, whereas the latters suppose fluctuation effects. According to  the cross-linking process is not sensible to strain $`\epsilon `$ and corresponding localization parameter $`\omega `$ is proportional to the cross-link number $`N`$ but not to the difference $`NN_c`$, as in Eqs.(25), (26). In our opinion, this is caused by non-self-consistency of the approach  in sense that the stress field $`\sigma `$ is switched off. On the other hand, making use of the statistical scheme  arrives at the strain-dependence for the localization parameter $`\omega `$ due to appearance of the polymer network entanglement, whereas within the above phenomenological approach this dependence has to be postulated.

A contribution to the elastic energy caused by the external strain is determined by equality

$$\mathrm{\Delta }E(e)|E\left(\epsilon (e)\right)E(\epsilon _0)|.$$
(30)
With use of Eqs. (22), (27), (29) and expansion $`\mathrm{ln}(1+x)xx^2/2+x^3/3`$$`x1`$ the expression

$$\mathrm{\Delta }E(e)N_0\frac{Tϵ^3}{2}e^2$$
(31)
is easily obtained where only the first non-zero term is kept. Comparing this relation with the usual expression for the elastic energy

$$\mathrm{\Delta }E(e)V\frac{G}{2}e^2$$
(32)
where $`V`$ is volume, we arrive at the final expression for the shear modulus of the gel state of the polymer network:

$$G=\frac{Tϵ^3}{\mathrm{\Omega }},ϵ\frac{N_0N_c}{N_c}1,\mathrm{\Omega }\frac{V}{N_0}.$$
(33)
The notable peculiarity of this result consists in that, in accordance with previous result (11), the shear modulus is proportional to the third power of the distance $`ϵ`$. It is worthwhile to note that an expression of this kind can be obtained only within framework of the synergetic approach used, but not on the basis of the phase transition theory by Landau.

Finally, comparing Eqs. (11) and (33), we obtain the relation between the micro- and macroscopic parameters of the gel under consideration in the proximity of the critical point ($`ϵ1`$):

$$G=g\frac{\xi }{𝒩}\frac{T}{\mathrm{\Omega }},g\frac{3}{2}f^1\frac{(f2)^2}{(f1)^4}.$$
(34)


format: text

text:
# Can a Very Low-Luminosity and Cold White Dwarf be a Self-gravitating Bose Condensed System

## Abstract

An entirely new model for the structure as well as for the cooling mechanism of white dwarfs has been proposed. We have argued that the massive part of the constituents of white dwarfs- the positively charged ions are boson and under the extreme physical condition (density and temperature) at the interior, it is possible to have condensation of charged bose gas. We have tried to establish that a cold white dwarf is a self gravitating charged bose condensed system.

PACS:97.20.Rp, 97.10.Cv, 03.75.F, 05.30.Jp

Unlike the other branches of physical science, astrophysics deals with giant macroscopic objects and also tiny constituents of microscopic world. In the astrophysical calculations very often we use the numbers as large as Planck mass ($`M_p`$) and very small number like Planck’s constant ($`h`$). The observational data from most of the giant macroscopic stellar objects, e.g., stars, white dwarfs, neutron stars, etc. are analyzed by various models based on micro-physical processes. To be more precise, the processes like stellar evolutions, supernova explosions, structure of white dwarfs and neutron stars etc. are explained by microscopic models based on nuclear physics, hydrodynamics, quantum mechanics, statistical mechanics, etc. Various astrophysical processes, e.g., emission of high energy radiations (in particular X-rays and $`\gamma `$-rays), neutrinos, many physical phenomena from pulsar observations, viz., pulsar glitches, existence of strong magnetic fields in neutron stars, etc. are also explained by various micro-physical processes. The theoretical analysis of thermal evolution and also the evolution of magnetic fields of neutron stars are based on the physics of superfluidity, superconductivity etc. of neutron matter . Similarly, the theoretical studies on the structure and the thermal evolution (which is again related to luminosity function) of white dwarfs of different masses and temperatures, in which we are especially interested in this article are also based on some interesting micro physical models . It is believed that the structure of white dwarfs is built on the following two assumptions. Electrons are all degenerate and they contribute to the total pressure of the system. On the other hand, ions are non-degenerate and almost the entire mass of the star comes from their rest masses. Based on these two assumptions and by employing the condition of hydrostatic equilibrium, the equation for the structure of white dwarfs can be obtained. Now the basic constituents of white dwarfs are mainly oxygen, carbon, helium with an envelope of hydrogen gas. It is however, sometime assumed that the whole star is made up of either carbon or helium or a composite structure of these three elements, with the heaviest at the core and lightest at the crust. The structure of a white dwarf is strictly determined by its mass. In some relatively massive white dwarfs, one can think of the presence of heavier elements like neon, magnesium or even iron within the stars. The class of white dwarfs which have originated as main sequence stars with masses around $`5M_{}`$ or below are carbon-oxygen ($`CO`$) white dwarfs, the other category, which originated from main sequence stars with masses from $`5M_{}`$ to about $`12M_{}`$ are oxygen-neon-magnesium ($`ONeMg`$) white dwarfs. It is also possible to have white dwarfs with helium rich interiors. Stars of low enough mass are supposed to evolve, after its hydrogen burning phase, to a stage where degeneracy appears before helium burning. But such isolated stellar objects are not expected to have yet evolved beyond the main sequence. However, it is possible to have white dwarfs with He-core  from more massive progenitor in a time scale smaller than the age of the galaxy if they belong to a binary system which allows mass loss before He-burning. Now because of high density all these elements within the white dwarfs are in completely ionized condition. The stability of the star against gravitational collapse is supported by degenerate electron pressure. Now the conventional model of white dwarf cooling is based on the following three consecutive processes: (i) relaxation of lattice thermal vibration, (ii) crystallization of the ionic part which are in the gaseous/liquid phase and (iii) further release of energy by the separation of different masses during crystallization in presence of gravity  (see also ). The total available energy in these processes depends on the constituent mass or equivalently on the mass of the white dwarf and temperature of the system. The time scale for cooling process is $`10`$Gyr . In this article we shall try to establish an entirely new picture for both the structure and the cooling of white dwarfs based on the condensation of charged bose gas. Now it is well known from the nuclear structure studies, that the kind of isotopes of He, C, O, Mg, Ne, etc. present in white dwarfs are all bosons. The white dwarf at relatively high temperature is therefore a self gravitating (one-component or multi-component) positively charged giant macroscopic bose system with a back ground of electron gas (negatively charged degenerate Fermi system). Based on this fact, we have tried to establish from a very simple physical argument that within the white dwarf it is possible to achieve a temperature which leads to condensation of charged bose gas. We shall also try to explain the cooling of white dwarfs with this new physical picture by assuming the condensation of bose gas as a first order phase transition.

Before we go into the detail discussion of bose condensation in a cold white dwarf, we give a very brief introduction on the condensation of ideal and weakly interacting hard-core bose gases. The condensation of weakly interacting bose gas- one of the oldest subjects has got a new dimension after the remarkable experimental achievement of bose condensation of alkali atoms in magnetic traps in the laboratory (, also see  for an overall knowledge in this subject). This discovery has given re-birth to this old subject and created an enormous interest in both theoretical and experimental studies of weakly interacting bose gas . It is well known that in the case of an ideal bose gas, the condition determining the fugacity $`z`$ in the normal phase is $`\lambda ^3n=g_{3/2}(z)`$, where $`n`$ is the particle density ,

$$\lambda =\left(\frac{2\pi \mathrm{}^2}{mkT}\right)^{1/2}$$
(1)
is the thermal wavelength and $`g_\nu (z)=_{l=0}^{\mathrm{}}z^\nu /l^\nu `$. Now the condition that the zero momentum state (i.e. the condensed state) is macroscopically occupied is $`\lambda ^3n>\zeta (3/2)`$, where $`\zeta (3/2)=g_{3/2}(z=1)2.612`$, the standard zeta function. Hence the transition temperature for the condensation of ideal bose gas is given by,

$$T_0=\frac{2\pi \mathrm{}^2}{mk}\left[\zeta \left(\frac{3}{2}\right)\right]^{\frac{2}{3}}n^{\frac{2}{3}}$$
(2)
On the other hand for an imperfect bose gas with non negligible hard core radius of the constituent, it was shown in several interesting theoretical calculations in the last few years, that the difference $`T=T_cT_00`$, where $`T_c`$ is the critical temperature in the case of imperfect bose gas. But there is no unique relation for the variation of $`T`$ with scattering length $`a`$. Not even there is consensus on the sign of $`T`$. Following the argument that a spatial repulsion of the constituents is equivalent to the momentum space attraction- which leads to condensation of hard core bose gas at relatively high temperature. Now for a system of identical bosons, assuming that the interaction potential acts locally (i.e. the range of interaction is small compared to the interparticle distance) and characterized entirely by $`s`$-wave scattering length $`a`$, the limit for which the quantum many body perturbation is valid is $`a\lambda `$. With these assumptions it was shown  by Stoof that $`Ta^{\frac{3}{2}}`$, Bijlsma and Stoof showed that $`Ta^{\frac{1}{2}}`$. The Monte-Carlo simulation gives $`T=T_0C_0(na^3)^\gamma `$. Where $`C_0=0.34\pm 0.06`$ and $`\gamma =0.34\pm 0.003`$. The mean field calculation shows $`T=0.7T_0(na^3)^{\frac{1}{3}}`$. The simple calculation by Huang shows that $`T=3.527T_0(na^3)^{\frac{1}{6}}`$. All these interesting results show that although there is no consensus on how $`T`$ depends on scattering length, however, there is no confusion that the condensation of imperfect bose gases occur at relatively higher temperature compared to that of an ideal bose gas.

We shall now go back to the discussion on the possibility of bose condensation in cold white dwarfs. The system we have considered contains massive ions carrying $`Z`$-unit positive charge ($`Z`$ is fixed for one-component systems, on the other hand it is a variable in the case of multicomponent systems) which form a non-degenerate charged bose gas neutralized by a back ground of degenerate electron gas. Now in gases with short range forces the weak coupling limit corresponds to low density range, whereas the weak coupling limit in the case of charged bose gases interacting via coulomb force, this corresponds to the high density limit. Hence in the case of white dwarfs of mass density $`10^6`$gm/cc, and at temperature $`10^610^8`$K, the weak coupling limit is a good approximation. Therefore, if we simulate the screened coulomb potential of a positive ion by some equivalent hard-core potential with range $`r_sa`$, the scattering length, we can apply the results of hard-core bose gas in this dense bosonic system. The scattering length $`a`$ defined above indicates the beginning of the repulsive region of the ion.

To get an idea of the order of magnitude of transition temperature in such dense bosonic systems we equate the thermal de Broglie wavelength (eqn.(1)) with the inter-ionic spacing $`n^{1/3}`$. The transition temperature $`T_c^{(1)}`$ in $`{}_{}{}^{o}K`$ obtained from this equality is displayed in Table I for three different elements (He, C and O). We have considered four possible mass densities. The transition temperatures so obtained for various physical conditions are found to be well within the limit of expected internal temperature of white dwarfs. Therefore the possibility of bose condensation inside cold white dwarfs can not be ruled out. The critical temperature $`T_c^{(0)}`$ obtained from eqn.(2), assuming ideal bose gas model is also shown in the table for the same set of elements and densities. We have also displayed the modified value of critical temperature $`T_c^{(h)}`$ from Huang’s paper  (the values do not change significantly with other results as cited before). Now it is almost impossible to obtain the scattering length $`a`$ within the white dwarf matter, we therefore treat it as a fixed parameter. Since we would like to use the results derived for hard-core bose gas obtained from perturbative calculation, keeping in mind the limit of its validity, we use the value $`a=0.01\lambda `$. We believe that even under extreme condition ($`\rho =10^9`$gm/cc), the scattering length can not go below this small value. The thermal wavelengths for two different cases: with and without hard-core radius are also presented in the table in $`\AA `$ unit. We have denoted these two quantities by $`\lambda _T^{(h)}`$ and $`\lambda _T^{(0)}`$ respectively. Assuming that the bose condensation is a first order phase transition, we have calculated the total amount of energy released in the process. The total available energy is given by $`E=LM`$, where $`LkT_c`$, the latent heat of condensation and $`M`$ is the total mass of the system. Assuming $`5.0\times 10^3`$Km as the radius of a typical white dwarf which undergoes a bose condensation transition at a certain critical temperature (for hard-core gas), we have shown the amount of energy liberated during this process within the star for the same set of elements and mass densities as discussed above. The detail analysis of thermal evolution of white dwarfs taking bose condensation into account will be presented in some future publication.

From the tabulated results it may be inferred that the critical temperature for condensation increases with density. The transition process will therefore be favored in massive white dwarfs. It is also obvious that for a given density the lighter ions condense at relatively higher temperature. Now the ionic distribution in a multicomponent white dwarf is such that the matter at the core region is mainly dominated by massive ions, whereas the lighter ions built the crustal part. Therefore, the argument given by Lamb  long ago on the possibility of separation of normal and bose condensed phases in presence of gravity may be ruled out. However, the process strongly depends on so many factors, e.g., the mass, central density and temperature, the basic constituents and finally on the radial distribution of mass and temperature within the star. If there is a phase separation in presence of gravity, some extra energy will also be released in the process. A detail analysis will be presented in a future paper. Therefore the final conclusion is that the possibility of bose condensation can not be ruled out, particularly in massive white dwarfs of very low luminosity, at least at the central region. The structure of such massive white dwarfs could therefore be a bose condensed core with crystalline normal crustal matter. We can compare such structures with the typical structure of a neutron star with superfluid core and normal neutron matter crust. In a future publication, we shall discuss the most important issue- the controversy of white dwarf ages in globular clusters using this new ideas of white dwarf structure and cooling mechanism introduced in this letter.

| Element | Helium | | | |
| --- | --- | --- | --- | --- |
| $`\rho `$ (gm/cc) | $`10^6`$ | $`10^7`$ | $`10^8`$ | $`10^9`$ |
| $`T_c^{(1)}`$oK | $`6.8\times 10^5`$ | $`3.1\times 10^6`$ | $`1.5\times 10^7`$ | $`6.8\times 10^7`$ |
| $`T_c^{(0)}`$oK | $`1.1\times 10^5`$ | $`5.3\times 10^5`$ | $`2.4\times 10^6`$ | $`1.1\times 10^7`$ |
| $`T_c^{(h)}`$oK | $`2.6\times 10^5`$ | $`1.2\times 10^6`$ | $`5.7\times 10^6`$ | $`2.6\times 10^7`$ |
| $`\lambda _T^{(0)}`$$`\AA `$ | $`2.6\times 10^2`$ | $`1.2\times 10^2`$ | $`5.6\times 10^3`$ | $`2.6\times 10^3`$ |
| $`\lambda _T^{(h)}`$$`\AA `$ | $`1.7\times 10^2`$ | $`7.9\times 10^3`$ | $`3.7\times 10^3`$ | $`1.7\times 10^3`$ |
| $`\mathrm{\Delta }E`$ (ergs) | $`1.0\times 10^{46}`$ | $`5.0\times 10^{47}`$ | $`2.0\times 10^{49}`$ | $`1.0\times 10^{51}`$ |
| Element | Carbon | | | |
| $`\rho `$ (gm/cc) | $`10^6`$ | $`10^7`$ | $`10^8`$ | $`10^9`$ |
| $`T_c^{(1)}`$oK | $`1.1\times 10^5`$ | $`5.0\times 10^5`$ | $`2.3\times 10^6`$ | $`1.1\times 10^7`$ |
| $`T_c^{(0)}`$oK | $`1.8\times 10^4`$ | $`8.5\times 10^4`$ | $`3.9\times 10^5`$ | $`1.8\times 10^6`$ |
| $`T_c^{(h)}`$oK | $`4.2\times 10^4`$ | $`2.0\times 10^5`$ | $`9.1\times 10^5`$ | $`4.2\times 10^6`$ |
| $`\lambda _T^{(0)}`$$`\AA `$ | $`3.7\times 10^2`$ | $`1.7\times 10^2`$ | $`8.0\times 10^3`$ | $`3.7\times 10^3`$ |
| $`\lambda _T^{(h)}`$$`\AA `$ | $`2.5\times 10^2`$ | $`1.1\times 10^2`$ | $`5.3\times 10^3`$ | $`2.5\times 10^3`$ |
| $`\mathrm{\Delta }E`$ (ergs) | $`6.0\times 10^{44}`$ | $`2.8\times 10^{46}`$ | $`1.7\times 10^{48}`$ | $`6.1\times 10^{49}`$ |
| Element | Oxygen | | | |
| $`\rho `$ (gm/cc) | $`10^6`$ | $`10^7`$ | $`10^8`$ | $`10^9`$ |
| $`T_c^{(1)}`$oK | $`6.7\times 10^4`$ | $`3.1\times 10^5`$ | $`1.4\times 10^6`$ | $`6.7\times 10^6`$ |
| $`T_c^{(0)}`$oK | $`1.1\times 10^4`$ | $`5.2\times 10^4`$ | $`2.4\times 10^5`$ | $`1.1\times 10^6`$ |
| $`T_c^{(h)}`$oK | $`2.6\times 10^4`$ | $`1.2\times 10^5`$ | $`5.6\times 10^5`$ | $`2.6\times 10^6`$ |
| $`\lambda _T^{(0)}`$$`\AA `$ | $`4.1\times 10^2`$ | $`1.9\times 10^2`$ | $`8.9\times 10^3`$ | $`4.1\times 10^3`$ |
| $`\lambda _T^{(h)}`$$`\AA `$ | $`2.7\times 10^2`$ | $`1.3\times 10^3`$ | $`5.8\times 10^3`$ | $`2.8\times 10^3`$ |
| $`\mathrm{\Delta }E`$ (ergs) | $`3.0\times 10^{44}`$ | $`1.3\times 10^{46}`$ | $`6.0\times 10^{47}`$ | $`2.8\times 10^{49}`$ |

TABLE I

ACKNOWLEDGEMENT: The work is supported by The Department of Science & Technology, Government of India, sanction no.: SP/S2/K3/97(PRU).


text:
# Modelling of dislocation generation and interaction during high-speed deformation of metals.Presented at the Third Workshop on High-speed Plastic Deformation, Hiroshima Institute of Technology, August 2000.

## 1 Introduction

In a recent paper, Kiritani *et al.* reported that a large number of vacancies were produced during high-speed heavy plastic deformation of thin foils of fcc metals . They observed a large density of stacking-fault tetrahedra but very low dislocation densities in the foils after deformation. As a possible explanation, they propose a dislocation-free deformation mechanism.

In this paper we investigate the deformation processes in metals under very high strain rates. Our approach is that of Molecular Dynamics (MD), where Newton’s second law is solved numerically for all the atoms in the system under consideration. The main advantage of this method is that no *a priori* assumptions are made about the deformation mechanism. Structures such as dislocations appear automatically from the motion of the atoms. The simulation thus becomes “unbiased”, i.e. the result is independent of our prior assumptions about the active deformation mechanisms.

We have chosen to simulate the deformation processes in systems that are initially dislocation-free (i.e. without conventional dislocation sources such as Frank-Read sources). We start with single crystals in order to avoid dislocation emission from grain boundaries. Simulations of high-speed deformation of polycrystalline (nanocrystalline) material have been reported elsewhere . Although many metals were studied experimentally by Kiritani et al., we have chosen to focus our simulations on one metal: copper.

## 2 Simulation methods and setup

In order to get meaningful results from a molecular dynamics simulation, it is essential that the interatomic forces used in the simulation are a good approximation to the forces the atoms would experience in a real physical system with the same atomic configuration. We use the effective medium theory (EMT)  to describe the interatomic interactions, as the EMT potential for copper has been well tested.

The simulated systems were single crystals. Two different system sizes were used, one with approximately 95 000 atoms and one with approximately 765 000 atoms. The linear dimensions of the two systems are approximately 10.5 nm and 21 nm. In both cases periodic boundary conditions were applied to the simulation, so the system effectively becomes an infinite crystal.

In the smaller system four vacancies were created in the initial configuration, by removing four random atoms. They were introduced in the unlikely case that vacancies were somehow actively participating in the deformation. The system is deformed along the  direction at a temperature of 450 K. We have chosen this temperature since we consider it likely that significant local heating is occurring in the experiments due to the high strain rate.

As the preexisting vacancies in the smaller system did not participate in the deformation, no vacancies were introduced in the larger system. It was deformed along the \[$`\overline{1}05`$\] direction at the same temperature.

The systems were deformed at a constant strain rate, while keeping the stress in the transverse direction approximately constant. This procedure is further described in ref. . In all cases a strain rate of $`10^9s^1`$ was used (compared to approximately $`10^5s^1`$ in the experiments). To keep the temperature of the system approximately constant during the simulation, Langevin dynamics  are used, i.e. a friction and a fluctuating force are added to the equations of motion of the atoms. We use a timestep of 5 fs, safely below the value where the dynamics becomes unstable.

During the simulations, configurations were rapidly quenched, and the local crystal structure was identified by common neighbor analysis . This was used to generate plots where all atoms except atoms at the dislocation cores were made invisible, allowing the visualization of the dislocation structures.

## 3 Results

Figure 1 shows the deformation of the smaller system. In figure 1(b) the first dislocation activity appears in the form of a few dislocation loops. The loops are nucleated homogeneously, i.e. without the presence of dislocation sources or crystal defects. The exact time and position of the nucleation event is determined by random thermal fluctuations. The loops are faulted loops consisting of a single partial dislocation with a stacking fault inside the loop. The loops expand rapidly. Due to the periodic boundary conditions the dislocations cannot disappear from the sample. They will continue to move causing significant plastic deformation. As the dislocations on different glide planes interact, dislocation reactions occur and a number of triple junctions appear. These typically consist of two Shockley partials that meet to form a stair-rod dislocation. Continuous formation and destruction of these sessile dislocations are seen in the simulations.

The dislocations observed in the simulations are mainly single Shockley partials. The homogeneous dislocation nucleation mechanism makes it very difficult to nucleate the second partial. Dislocations are only nucleated at extremely high stresses, but as soon as the first partial dislocation is created, it locally shields the stress field and prevents the second partial from nucleating. It probably also plays a role that the stacking fault energy is low.111Copper has a relatively low stacking fault energy, and furthermore the stacking fault energy in the simulation is approximately a factor of two below the experimental value (30 mJ/m2 versus 52 mJ/m2). At high stresses the energy of the stacking faults become small compared to other energies in the system, such as the energy released by a moving dislocation.

After some dislocation activity has occurred we observe the first creation of vacancies in figure 1(c). A “sausage-like” string of vacancies has appeared in the middle of the system. Further dislocation activity breaks up the string of vacancies into smaller clusters.

To analyse the process by which such a cluster of vacancies is generated, we repeated the relevant part of the simulation, using a higher time resolution for the analysis of the process. This revealed the detailed process generating these vacancies. The mechanism is illustrated in figure 2. Dislocation interactions have created dislocations where segments are on different glide planes. This makes it possible for a segment of a dislocation to be annihilated. If two edge dislocations gliding on adjacent glide planes annihilate, a string of vacancies is formed.222If the dislocations had the opposite signs, a string of interstitials should be formed. This is not observed in the simulations, probably because the energy of such a string is too high, so the dislocations pass each other instead of annihilating. This appears to be the main source of vacancies in the simulations.

At the end of the simulation (figure 1(d)) a large number of vacancies have been formed, apparently by this mechanism.

The larger system behaves in a similar fashion. Figure 3 shows the larger system after the deformation. As a large number of dislocations were present, it would have been very difficult to see anything in the figure if the entire system were plotted. For this reason only a slice through the system is shown. Again a large number of both dislocations and vacancies are seen.

## 4 Discussion

The simulations presented here indicate that the deformation of fcc metals at ultra-high strain rate occurs by mechanisms that are not very different from what is seen at more moderate strain rates. Even in a situation where there are no preexisting dislocations, and no dislocation sources, dislocations are nucleated and are the main carriers of the deformation. The main difference between the deformation mechanism in the simulations and in slowly deforming materials is how the dislocations are nucleated. This difference may not be present in a real physical sample under ultra-high-speed deformation, since dislocation sources will be present, even though it is possible that other types of dislocation sources are activated under high strain rate conditions.

The high concentration of vacancies in the simulations is due to the very high dislocation density during the high-speed deformation. As the vacancies are generated when opposite-signed dislocations collide, it must be expected that the generation rate depends on the dislocation densities obtained and not just on the total deformation. It is therefore expected that much higher vacancy densities should be observed experimentally after ultra-rapid deformation than after deformation with moderate strain rate.

In the simulations presented here there are no free surfaces, and no efficient way of removing dislocations except by annihilation. We therefore expect a higher dislocation density and a higher vacancy production rate than those which would be observed with a more realistic geometry.

To test the effect of free surfaces, we have made a similar simulation with four free surfaces (i.e. the system has the form of a rod being pulled). In this case dislocations were nucleated at the corners, passed through the system and disappeared at the free surfaces on the other side of the sample. Only a few dislocations were moving through the system at the same time. No vacancies were generated, as the dislocations never collided. Preliminary simulations with two free surfaces (i.e. foil geometry) indicate that this geometry behaves much like the bulk simulations, as dislocations perpendicular to the surfaces cannot move to the surface and disappear. A large dislocation density builds up, and vacancies are generated.

## 5 Conclusions

Atomic-scale computer simulations using molecular dynamics indicate that ultra-high-speed deformation of fcc copper occurs by the nucleation and motion of dislocations. Even in the absence of dislocation sources, dislocations will be nucleated at sufficiently high stresses.

During the deformation, a large number of vacancies are formed, in agreement with experimental investigations of ultra-high strain rate deformation of thin foils of fcc metals . The vacancies are produced when segments of edge dislocations moving on adjacent slip planes annihilate. This generates “strings” of dislocations. Further dislocation activity tends to break these strings into smaller clusters. Finally, over time spans much longer than those which we can simulate using molecular dynamics, the vacancies will diffuse and will eventually cluster, forming stacking fault tetrahedra.

## 6 Acknowledgments

We would like to thank Prof. Kiritani for fruitful discussions. Center for Atomic-scale Materials Physics (CAMP) is sponsored by the Danish National Research Foundation. This work was done as a collaboration between CAMP and the Engineering Science Center for Structural Characterization and Modelling of Materials at Materials Research Department, Risø. Parallel computer time was financed by the Danish Research Councils through grant no. 9501775.


text:
# References

Non-Fermi liquid angle resolved photoemission lineshapes of Li0.9Mo6O17 A recent Letter by Xue *et al.*  reports observations from angle resolved photoemission spectroscopy (ARPES) on quasi one-dimensional (q-1D) Li0.9Mo6O17 that above its $`T_X`$ $``$$`24`$K transition a peak dispersing to define its Fermi surface (FS) develops Fermi energy ($`E_F`$) weight requiring a Fermi liquid (FL) description. This finding contradicts our report  of a non-FL lineshape in this material. The reasoning in  that this new finding was enabled by improved angle resolution is flawed. Rather, we point out here that the data of  have fundamental differences from other ARPES results and also band theory. Therefore, until the origin of these differences are learned, the claims of  must be held in abeyance. These claims also include the report of an 80meV gap below $`T_X`$, which contradicts the zero gap found in optical spectroscopy  and magnetic susceptibility and exceeds by more than two orders of magnitude the value (0.3meV) implied by a gap-model interpretation of the resistivity rise below $`T_X`$ .

Improved angle resolution is not relevant for the claimed FL lineshape. Because the k=kF lineshape for both the FL and the Tomonaga-Luttinger (TL) model (with $`\alpha <1`$) are singular at $`E_F`$, the $`E_F`$ weight for both increases steadily as the angle resolution improves. Indeed, it is well known that one must angle-integrate (k-sum) to test for the surprising difference between the Fermi edge of the FL and the $`E_F`$ power law onset of the TL model. Xue *et al*. sum ARPES data along the q-1D $`\mathrm{\Gamma }`$–Y direction over $`\mathrm{\Delta }𝐤`$=0.2Å-1 and report a Fermi edge, whereas we always find only a power law onset at $`E_F`$ in angle summed spectra, including our new high resolution spectra shown below. This difference arises because the data are fundamentally different.

Fig. 1 shows various $`\mathrm{\Gamma }`$–Y data sets. The kF values show small variations among the data sets and so we place importance on the k value relative to kF. Panel (a) shows previously unpublished data taken by us at photon energy h$`\nu `$=24 eV on literally the same cleaved surface as for the data of . Panel (b) shows earlier data with h$`\nu `$=21.2 eV by Grioni *et al*.  The data sets for the two h$`\nu `$ values are consistent, both showing bands A through D in good basic agreement with band theory, as labeled. Only bands C,D cross $`E_F`$, becoming degenerate before the crossing. In addition, there is an unexplained tendency, seen also in other Mo bronzes, for non-dispersive weight to cling to the bottom of the band ($``$$`0.6`$eV). Panels (a) and (b) establish the consistency of the data of  and . The special k-path parallel to $`\mathrm{\Gamma }`$–Y in  was chosen because both bands C,D are especially strong all the way to $`E_F`$. The bands C,D in (a)–(c) along $`\mathrm{\Gamma }`$–Y are in basic agreement with those in  (short dashes in (a)), and, importantly, show non-FL lineshapes as do the data of . A difference occuring for the filled band B along the two paths is implied by band theory but is not important here.

Panel (c) shows h$`\nu `$=24 eV data taken at the Wisconsin Synchrotron Radiation Center PGM beamline with an SES-200 Scienta analyzer over a narrow k-range near $`𝐤_F`$ where D is already too weak to see. These data have angle and energy resolution comparable to that of Xue *et al.* and fully agree with the data of (a) and (b) apart from generally increased sharpness, and increased $`E_F`$ weight for $`𝐤_F`$. It is then meaningful to compare the data of (c) directly to the h$`\nu `$=21.4 eV data of Xue *et al.* , shown in panel (d). Relative to $`𝐤_F`$, the two k-ranges are nearly the same. It is obvious by inspection that the two data sets are globally different, for example by the absence in (d) of peaks A,B and by the presence in (d) of a non-dispersing feature E which interferes with the presumed C,D lineshapes, and has no counterpart in the other data or in band theory. This work was supported by the U.S. DoE at U. Mich. and Iowa State, and by the U.S. NSF at U. Mich. and U. Wisc.

G.-H. Gweon1, J.D. Denlinger1, J.W. Allen1, C.G. Olson2, H. Höchst3, J. Marcus4 and C. Schlenker4

*1Randall Lab. of Physics, U. of Michigan*, *Ann Arbor, MI 48109*. *2Ames Lab., Iowa State U., Ames, Iowa 50011. 3U. of Wisconsin*, *SRC, Stoughton, WI 53589. 4LEPES, CNRS, BP 166, F-38042 Grenoble Cedex 9, France.*


text:
# Evidence for the exponential distribution of income in the USA

## 1 Introduction

The study of income distribution has a long history. Pareto Pareto  proposed in 1897 that income distribution obeys a universal power law valid for all times and countries. Subsequent studies have often disputed this conjecture. In 1935, Shirras Shirras  concluded: “There is indeed no Pareto Law. It is time it should be entirely discarded in studies on distribution”. Mandelbrot Mandelbrot  proposed a “weak Pareto law” applicable only asymptotically to the high incomes. In such a form, Pareto’s proposal is useless for describing the great majority of the population.

Many other distributions of income were proposed: Levy, log-normal, Champernowne, Gamma, and two other forms by Pareto himself (see a systematic survey in the World Bank research publication Kakwani ). Theoretical justifications for these proposals form two schools: socio-economic and statistical. The former appeals to economic, political, and demographic factors to explain the distribution of income (e. g. Levy ), whereas the latter invokes stochastic processes. Gibrat Gibrat  proposed in 1931 that income is governed by a multiplicative random process, which results in a log-normal distribution (see also Montroll ). However, Kalecki Kalecki  pointed out that the width of this distribution is not stationary, but increases in time. Levy and Solomon Solomon  proposed a cut-off at lower incomes, which stabilizes the distribution to a power law.

In this paper, we propose that the distribution of individual income is given by an exponential function. This conjecture is inspired by our previous work DY , where we argued that the probability distribution of money in a closed system of agents is given by the exponential Boltzmann-Gibbs function, in analogy with the distribution of energy in statistical physics. In Sec. 2, we compare our proposal with the census and tax data for individual income in USA. In Sec. 3, we derive the distribution function of income for families with two earners and compare it with the census data. The good agreement we found is discussed in Sec. 4. Speculations on the possible origins of the exponential distribution of income are given in Sec. 5.

## 2 Distribution of individual income

We denote income by the letter $`r`$ (for “revenue”). The probability distribution function of income, $`P(r)`$, (called the probability density in book Kakwani ) is defined so that the fraction of individuals with income between $`r`$ and $`r+dr`$ is $`P(r)dr`$. This function is normalized to unity (100%): $`_0^{\mathrm{}}P(r)𝑑r=1`$. We propose that the probability distribution of individual income is exponential:

$$P_1(r)=\mathrm{exp}(r/R)/R,$$
(1)
where the subscript 1 indicates individuals. Function (1) contains one parameter $`R`$, equal to the average income: $`_0^{\mathrm{}}rP_1(r)𝑑r=R`$, and analogous to temperature in the Boltzmann-Gibbs distribution DY .

From the Survey of Income and Program Participation (SIPP) census , we downloaded the variable TPTOINC (total income of a person for a month) for the first “wave” (a four-month period) in 1996. Then we eliminated the entries with zero income, grouped the remaining entries into bins of the size 10/3 k$, counted the numbers of entries inside each bin, and normalized to the total number of entries. The results are shown as the histogram in Fig. 1, where the horizontal scale has been multiplied by 12 to convert monthly income to an annual figure. The solid line represents a fit to the exponential function (1). In the inset, plot A shows the same data with the logarithmic vertical scale. The data fall onto a straight line, whose slope gives the parameter $`R`$ in Eq. (1). The exponential law is also often written with the bases 2 and 10: $`P_1(r)2^{r/R_2}10^{r/R_{10}}`$. The parameters $`R`$, $`R_2`$ and $`R_{10}`$ are given in line (c) of Table 1.

Plot B in the inset of Fig. 1 shows the data from the Panel Study of Income Dynamics (PSID) conducted by the Institute for Social Research of the University of Michigan Michigan . We downloaded the variable V30821 “Total 1992 labor income” for individuals from the Final Release 1993 and processed the data in a similar manner. Shown is the cumulative probability distribution of income $`N(r)`$ (called the probability distribution in book Kakwani ). It is defined as $`N(r)=_r^{\mathrm{}}P(r^{})𝑑r^{}`$ and gives the fraction of individuals with income greater than $`r`$. For the exponential distribution (1), the cumulative distribution is also exponential: $`N_1(r)=_r^{\mathrm{}}P_1(r^{})𝑑r^{}=\mathrm{exp}(r/R)`$. Thus, $`R_2`$ is the median income; 10% of population have income greater than $`R_{10}`$ and only 1% greater than $`2R_{10}`$. The points in the inset fall onto a straight line in the logarithmic scale. The slope is given in line (a) of Table 1.

The points in Fig. 2 show the cumulative distribution of tax returns vs income in 1997 from column 1 of Table 1.1 of Ref. Pub1304 . (We merged 1 k$ bins into 5 k$ bins in the interval 1–20 k$.) The solid line is a fit to the exponential law. Plot A in the inset of Fig. 2 shows the same data with the logarithmic vertical scale. The slope is given in line (e) of Table 1. Plot B in the inset of Fig. 2 shows the distribution of individual income from tax returns in 1993 SailerWeber . The logarithmic slope is given in line (b) of Table 1.

While Figs. 1 and 2 clearly demostrate the fit of income distribution to the exponential form, they have the following drawback. Their horizontal axes extend to $`+\mathrm{}`$, so the high-income data are left outside of the plots. The standard way to represent the full range of data is the so-called Lorenz curve (for an introduction to the Lorenz curve and Gini coefficient, see book Kakwani ). The horizontal axis of the Lorenz curve, $`x(r)`$, represents the cumulative fraction of population with income below $`r`$, and the vertical axis $`y(r)`$ represents the fraction of income this population accounts for:

$$x(r)=_0^rP(r^{})𝑑r^{},y(r)=\frac{_0^rr^{}P(r^{})𝑑r^{}}{_0^{\mathrm{}}r^{}P(r^{})𝑑r^{}}.$$
(2)
As $`r`$ changes from 0 to $`\mathrm{}`$, $`x`$ and $`y`$ change from 0 to 1, and Eq. (2) parametrically defines a curve in the $`(x,y)`$-space.

Substituting Eq. (1) into Eq. (2), we find

$$x(\stackrel{~}{r})=1\mathrm{exp}(\stackrel{~}{r}),y(\stackrel{~}{r})=x(\stackrel{~}{r})\stackrel{~}{r}\mathrm{exp}(\stackrel{~}{r}),$$
(3)
where $`\stackrel{~}{r}=r/R`$. Excluding $`\stackrel{~}{r}`$, we find the explicit form of the Lorenz curve for the exponential distribution:

$$y=x+(1x)\mathrm{ln}(1x).$$
(4)
$`R`$ drops out, so Eq. (4) has no fitting parameters.

The function (4) is shown as the solid curve in Fig. 3. The straight diagonal line represents the Lorenz curve in the case where all population has equal income. Inequality of income distribution is measured by the Gini coefficient $`G`$, the ratio of the area between the diagonal and the Lorenz curve to the area of the triangle beneath the diagonal: $`G=2_0^1(xy)𝑑x`$. The Gini coefficient is confined between 0 (no inequality) and 1 (extreme inequality). By substituting Eq. (4) into the integral, we find the Gini coefficient for the exponential distribution: $`G_1=1/2`$.

The points in Fig. 3 represent the tax data during 1979–1997 from Ref. Petska . With the progress of time, the Lorenz points shifted downward and the Gini coefficient increased from 0.47 to 0.56, which indicates increasing inequality during this period. However, overall the Gini coefficient is close to the value 0.5 calculated for the exponential distribution, as shown in the inset of Fig. 3.

## 3 Income distribution for two-earners families

Now let us discuss the distribution of income for families with two earners. The family income $`r`$ is the sum of two individual incomes: $`r=r_1+r_2`$. Thus, the probability distribution of the family income is given by the convolution of the individual probability distributions Feller . If the latter are given by the exponential function (1), the two-earners probability distribution function $`P_2(r)`$ is

$$P_2(r)=_0^rP_1(r^{})P_1(rr^{})𝑑r^{}=\frac{r}{R^2}e^{r/R}.$$
(5)
The function $`P_2(r)`$ (5) differs from the function $`P_1(r)`$ (1) by the prefactor $`r/R`$, which reflects the phase space available to compose a given total income out of two individual ones. It is shown as the solid curve in Fig. 4. Unlike $`P_1(r)`$, which has a maximum at zero income, $`P_2(r)`$ has a maximum at $`r=R`$ and looks qualitatively similar to the family income distribution curves in literature Levy .

From the same 1996 SIPP that we used in Sec. 2 census , we downloaded the variable TFTOTINC (the total family income for a month), which we then multiplied by 12 to get annual income. Using the number of family members (the variable EFNP) and the number of children under 18 (the variable RFNKIDS), we selected the families with two adults. Their distribution of family income is shown by the histogram in Fig. 4. The fit to the function (5), shown by the solid line, gives the parameter $`R`$ listed in line (d) of Table 1. The families with two adults and more than two adults constitute 44% and 11% of all families in the studied set of data. The remaining 45% are the families with one adult. Assuming that these two classes of families have two and one earners, we expect the income distribution for all families to be given by the superposition of Eqs. (1) and (5): $`0.45P_1(r)+0.55P_2(r)`$. It is shown by the solid line in the inset of Fig. 4 (with $`R`$ from line (d) of Table 1) with the all families data histogram.

By substituting Eq. (5) into Eq. (2), we calculate the Lorenz curve for two-earners families:

$$x(\stackrel{~}{r})=1(1+\stackrel{~}{r})e^{\stackrel{~}{r}},y(\stackrel{~}{r})=x(\stackrel{~}{r})\stackrel{~}{r}^2e^{\stackrel{~}{r}}/2.$$
(6)
It is shown by the solid curve in Fig. 5. Given that $`xy=\stackrel{~}{r}^2\mathrm{exp}(\stackrel{~}{r})/2`$ and $`dx=\stackrel{~}{r}\mathrm{exp}(\stackrel{~}{r})d\stackrel{~}{r}`$, the Gini coefficient for two-earners families is: $`G_2=2_0^1(xy)𝑑x=_0^{\mathrm{}}\stackrel{~}{r}^3\mathrm{exp}(2\stackrel{~}{r})𝑑\stackrel{~}{r}=3/8=0.375`$. The points in Fig. 5 show the Lorenz data and Gini coefficient for family income during 1947–1994 from Table 1 of Ref. history . The Gini coefficient is very close to the calculated value 0.375.

## 4 Discussion

Figs. 1 and 2 demonstrate that the exponential law (1) fits the individual income distribution very well. The Lorenz data for the individual income follow Eq. (4) without fitting parameters, and the Gini coefficient is close to the calculated value 0.5 (Fig. 3). The distributions of the individual and family income differ qualitatively. The former monotonically increases toward the low end and has a maximum at zero income (Fig. 1). The latter, typically being a sum of two individual incomes, has a maximum at a finite income and vanishes at zero (Fig. 4). Thus, the inequality of the family income distribution is smaller. The Lorenz data for families follow the different Eq. (6), again without fitting parameters, and the Gini coefficient is close to the smaller calculated value 0.375 (Fig. 5). Despite different definitions of income by different agencies, the parameters extracted from the fits (Table 1) are consistent, except for line (e).

The qualitative difference between the individual and family income distributions was emphasized in Ref. SailerWeber , which split up joint tax returns of families into individual incomes and combined separately filed tax returns of married couples into family incomes. However, Refs. Pub1304  and Petska  counted only “individual tax returns”, which also include joint tax returns. Since only a fraction of families file jointly, we assume that the latter contribution is small enough not to distort the tax returns distribution from the individual income distribution significantly. Similarly, the definition of a family for the data shown in the inset of Fig. 4 includes single adults and one-adult families with children, which constitute 35% and 10% of all families. The former category is excluded from the definition of a family for the data history  shown in Fig. 5, but the latter is included. Because the latter contribution is relatively small, we expect the family data in Fig. 5 to approximately represent the two-earners distribution (5). Technically, even for the families with two (or more) adults shown in Fig. 4, we do not know the exact number of earners.

With all these complications, one should not expect perfect accuracy for our fits. There are deviations around zero income in Figs. 1, 2, and 4. The fits could be improved there by multiplying the exponential function by a polynomial. However, the data may not be accurate at the low end because of underreporting. For example, filing a tax return is not required for incomes below a certain threshold, which ranged in 1999 from $2,750 to $14,400 Pub1040 . As the Lorenz curves in Figs. 3 and 5 show, there are also deviations at the high end, possibly where Pareto’s power law is supposed to work. Nevertheless, the exponential law gives an overall good description of income distribution for the great majority of the population.

## 5 Possible origins of exponential distribution

The exponential Boltzmann-Gibbs distribution naturally applies to the quantities that obey a conservation law, such as energy or money DY . However, there is no fundamental reason why the sum of incomes (unlike the sum of money) must be conserved. Indeed, income is a term in the time derivative of one’s money balance (the other term is spending). Maybe incomes obey an approximate conservation law, or somehow the distribution of income is simply proportional to the distribution of money, which is exponential DY .

Another explanation involves hierarchy. Groups of people have leaders, which have leaders of a higher order, and so on. The number of people decreases geometrically (exponentially) with the hierarchical level. If individual income increases linearly with the hierarchical level, then the income distribution is exponential. However, if income increases multiplicatively, then the distribution follows a power law hierarchy . For moderate incomes below $100,000, the linear increase may be more realistic. A similar scenario is the Bernoulli trials Feller , where individuals have a constant probability of increasing their income by a fixed amount.

###### Acknowledgements.

We are grateful to D. Jordan, M. Weber, and T. Petska for sending us the data from Refs. Pub1304 , SailerWeber , and Petska , to T. Cranshaw for discussion of income distribution in Britain, and to M. Gubrud for proofreading of the manuscript.