Statistics of Binary Exchange of Energy or Money

Why does the Maxwell-Boltzmann energy distribution for an ideal classical gas have an exponentially thin tail at high energies, while the Kaniadakis distribution for a relativistic gas has a power-law fat tail? We argue that a crucial role is played by the kinematics of the binary collisions. In the classical case the probability of an energy exchange far from the average (i.e., close to 0% or 100%) is quite large, while in the extreme relativistic case it is small. We compare these properties with the concept of “saving propensity”, employed in econophysics to define the fraction of their money that individuals put at stake in economic interactions.


Introduction
In this work, we study with statistical methods a problem concerning the elastic scattering of two particles in the relativistic regime and we discuss the similarities between this physical process and the economic interactions occurring in the so-called wealth exchange models [1,2].A unifying element in this analysis is the generalized κ-distribution [3][4][5][6], which has proven to be helpful for the description of both kinds of phenomena.We start from the familiar assumption of statistical mechanics according to which the macroscopic distributions (of energy or income) are determined by the microscopic interactions (between molecules or, respectively, individuals).Such interactions can be of many different sorts, but we are mostly interested into properties which are independent from the details of the specific system; namely we are interested into the conditions under which a "fat" power-law tail can appear in the macroscopic distribution, instead of the thin exponential tail typical of the Boltzmann distribution.
In our computation we shall suppose that all the permitted final states of a binary collision (those compatible with conservation laws) have the same probability.For this reason we mention in the title the "statistics of binary exchange", and not the dynamics of this exchange."Statistics" is meant here in the same sense as in the term "Statistical Mechanics".In fact, it is known from experimental data on relativistic gases and cosmic rays that the fat tails in their energy distributions have some universal features; it is also known that Pareto tails in income distributions have existed in various epochs and in countries with very different economic systems.(For some very recent data and related analysis see [7].)This suggests that the formation mechanism of the tails is essentially kinematical.In [3] G. Kaniadakis has shown that the deformation of the distribution function introduced by the parameter κ emerges naturally within Einstein's special relativity, so that one can see the κ-deformation as a pure relativistic effect.To this end, one first proves that the κ-deformed sum of the momenta of two particles is the additivity law for the relativistic momenta (the κ-sum is used in a generalized entropy minimization procedure leading to the κ-distribution).By considering, in the framework of special relativity, the κ-statistics of an ensemble of identical particles, one therefore arrives at a distribution with a power-law tail without any assumption on dynamics, but using only kinematics.
From the physical point of view, what are the main differences between a "classical" collision and a strongly relativistic collision?The answer is not unique.The techniques of relativistic kinetic theory [8] allow us to write the Boltzmann equation in relativistic form and to conduct a similar analysis as done in the classical theory.Nevertheless, the approach by Kaniadakis et al. shows that there is still room for a better understanding of the basic concepts.In our opinion, it is important to consider also quantities that are not Lorentz-invariant, and are therefore often disregarded in the traditional relativistic kinetic theory.We will focus our attention on the non-invariant fraction R of energy which is transferred from one particle to the other during a binary collision between identical particles.If one looks at the collision in the center of mass system, R is trivially zero, because the conservation laws of the total energy and momentum impose that the energy is equally distributed between the two particles before and after the scattering.What changes in the collision, as seen from the center of mass system, is only the direction of the momenta.Any other change is entirely due to the choice of a different reference system.
Nevertheless, for a gas contained in a finite volume there exists a privileged reference system, namely the system where the gas is, on average, at rest.If we look at a scattering process as the analogue of an economic interaction, the fraction R of money exchanged by an individual is obviously an important quantity.We point out that in any realistic income distribution the majority of the individuals have low or middle income, therefore it is important to evaluate R for an interaction between a poor individual (slow molecule) and a rich one (fast molecule).It has been suggested by some authors [9] that Pareto tails are formed when the rich part of a population enacts an effective policy for preserving its wealth in the interactions.This attitude can also be described as a high "saving propensity".In [10] we have argued that something similar happens in relativistic collisions, but the technique employed did not allow us to carry the analogy very far.In this work we improve the technique for the calculation of the R factor and we obtain further results.In Section 2 we evaluate the probability distribution of energy transfer in an elastic collision.We solve numerically the conservation law for energy and compute the integral over final states parametrized as a one-parameter line in phase space.In Section 3 we give the results of agent-based simulations based on money exchange rules which are the analogues of the energy balance in the collisions.Section 4 comprises our conclusions.The Appendices contain some background material on relativistic energy and momentum (Appendix A) and wealth exchange models in econophysics (Appendix B) that can be helpful for an interdisciplinary audience.
The question raised at the beginning of the abstract of this paper is therefore answered in two steps.In the first step (Section 2) we show that the probability of final states of a binary collision is different in the case of a low-energy and high-energy collision.In the second step (Section 3) we pass "from micro to macro" by simulating a large number of collisions occurring with this probability, and we look at what kind of equilibrium distribution they lead to in a large population of molecules.

General Expression
In this Section we consider the relativistic scattering of two identical particles of mass m and compute the probability distribution of the fraction R of kinetic energy transferred in the collision, supposing that all final states are equally probable.The incoming particle is labelled as "Particle 1" and the other, which is initially at rest, as "Particle 2".The initial velocity of Particle 1 is directed along x and is denoted simply as v.The reference system can be chosen in such a way that the final momenta of the particles have components only along x and y, denoted as p 1x , p 1y , p 2x , p 2y .Momentum conservation in the directions x and y imposes that p 2y = −p 1y ; (1) where p is the initial momentum of Particle 1: In the following we use units such that c = 1.Energy conservation requires that where E is the total energy E = m 2 + p 2 and E 1 and E 2 are the energies of the particles after the collision.The collision is supposed to be elastic.The fraction R of the initial kinetic energy of Particle 1 which is transferred to Particle 2 after the collision is The conservation laws ( 1)-( 3), taken together, make up a system with 4 unknowns and 3 equations, which has infinite solutions dependent on one parameter, say t.It is convenient to choose as parameter p 1x , the final momentum of Particle 1 along x.The value of this momentum spans a finite interval, because the minimum it can take is 0 (when the full momentum of the incoming particle is transferred to Particle 2, like in a frontal billiard scattering), and the maximum is p (when there is no momentum transfer at all).In the four-dimensional Euclidean vector space of the final momenta of the two particles (p 1x , p 1y , p 2x , p 2y ), the one-parameter subset of the solutions of Equations ( 1)-( 3) can be written as where t ∈ [0, p] and g(t) is the solution of Equation ( 3), which is now rewritten as The subset ( 5) is a path whose invariant length is given by In order to evaluate the average value in the final state of any quantity F depending on p 1x , p 1y , p 2x , p 2y , according to the ensemble averaging procedure familiar from statistical mechanics, we must compute the weighted integral where P is the appropriate probability density distribution.In this case, since we suppose that all possible final states of the scattering process are equiprobable, the function P is a constant function with support along the path (5).Therefore the four-dimensional integrals reduce to line integrals and we obtain where F is evaluated, through the momenta, as a function of t: In particular, we shall take as F the energy ratio R defined in (4).Note that L and F do not depend on the choice of the parameter t.The factor which multiplies F in the integral ( 9) is the parameter-invariant infinitesimal length ds of the path.This implements the physical intuition that any infinitesimal invariant interval of possible final states with the same length ds has the same probability.
Note that Equation ( 4) can be seen as interpolating between the purely non-relativistic case for which R is given by p 2  2 /p 2 , and the purely massless relativistic case for which R depends on the square root of the non-relativistic ratio.This transition in the effective functional dependence of R on p 2  2 and p 2 (as compared with the classical case) reflects the role of mass as a form of "savings" in the relativistic mass shell constraint.(We thank an anonymous referee for this remark.)

Numerical Solution of the Conservation Laws
Our strategy is to use the conservation laws ( 1)-( 3) to express the final momenta p 1y , p 2x , p 2y as functions of p 1x .So p 1x becomes our integration parameter t.However, the energy conservation relation (3) cannot be solved analytically.We thus transform it into a first-order differential equation by taking the total differential of Equation ( 3), which gives From this, with some straightforward algebraic steps, we obtain With the notation introduced in Equation ( 5) this can be rewritten in a form which is mathematically more transparent, recalling that m and p are fixed constants: In the following, however, we shall denote the parameter t by p 1x , according to its physical meaning.
We are interested into the numerical solution of this equation via a Runge-Kutta algorithm in the interval p 1x ∈ [0, p].It can be shown that the right hand sidediverges when p 1x → 0 + or p 1x → p − , so it is impossible to use the points p 1x = 0 or p 1x = p to set the the initial conditions of the numerical integration.It is possible, however, to impose the initial conditions at the value of p 1x which solves Equation (3) in the special case p 1x = p 1y ≡ ξ.In this case Equation (3) or its equivalent (6) become and can be solved numerically, thus finding appropriate initial conditions for the differential equation.
The solution of the latter can then be extended on the whole interval [0, p] by considering that: (1) The solution is symmetric with respect to p/2.(2) The solution in the restricted interval [ε, p − ε] converges when ε → 0. (One obtains numerical evidence of this property using small values of ε; the convergence follows from the fact that for any fixed ε the derivative g behaves as 1/ε α , with 0 < α < 1, when p 1x = ε or p 1x = p − ε.)An example of the solution, which we shall denote as g(p 1x ), is given in Figure 1.12), for a strongly relativistic case (p/m = γv = 100 γ; v 1).We use units c = 1 throughout the paper, and in this plot we also take m = 1 for graphical reasons.

Monte Carlo Evaluation of the Integral on Final States
The line integral (9) can be now computed through a Monte Carlo technique, since we have a numerical relation between p 1x and p 1y , while p 2x and p 2y can be obtained through the momentum conservation laws (1) and (2).For any function F of the final momenta, the average value F will only depend on the initial momentum p.We can therefore study the dependence of F on p distinguishing between the classical limit (γ → 1 + , v 1 in our units) and the strongly relativistic limit (γ 1, v → 1 − ).If we take F = R, it turns out that R is actually independent from p and is always equal to 1/2.This is due to the symmetry of the function g(p 1x ) (compare Figure 1).It is more interesting to consider the probability distribution of R in the interval I 1 = {p 1x ∈ [p/2, p]}.This corresponds to collisions in which Particle 2 has obtained from Particle 1 a kinetic energy between 0 and 1/2 of its initial kinetic energy (the values 0 and 1/2 being respectively obtained when p 1x = p and p 1x = p/2).The complementary interval in the possible final states, namely I 2 = {p 1x ∈ [0, p/2]}, represents the complementary cases of energy re-distribution, but since the two particles are identical, those cases are not distinguishable from the corresponding cases in I 1 , as far as their effect on the income distribution is concerned.
Let us make a numerical example with reference to Figure 2, which represents the probability distribution of R in a collision with low energy; then we shall compare this to the distribution in the extreme-relativistic case.The distribution in Figure 2 is obtained as follows: 1.
A numerical value of p/m is chosen, in this example p/m = 0.01, and Equation ( 12) is solved numerically, obtaining p 1y as a function g of p 1x : p 1y = g(p 1x ).

2.
A large number of random values of p 1x (typically 10 6 or 10 7 ) is generated in the interval I 1 .

3.
The corresponding value of R is calculated from the definition (4) using the function g(p 1x ) and Equations ( 1) and (2).

5.
In the graph, which is actually a histogram, the vertical axis represents, for 50 discretized values of R in the interval 0 ≤ R ≤ 1/2, the weighted number of counts, i.e., the number of times R has taken that value, multiplied by the weight.The part of the histogram for 1/2 ≤ R < 1 is symmetrical, and not shown.
For instance, in the almost-classical case of Figure 2 (p/m = 0.01), the probability that R = 0.4 = 40%, i.e., Particle 1 has lost to Particle 2 the 40% of its initial kinetic energy, is approximately 1/3 of the probability that R = 0.05 = 5%.Note that the probability of R = 60% is equal to the probability of R = 40% and in both cases the collision will have the same effective result, with one of the particles having 60% of the initial energy and the other 40%.
On the other hand, in the extreme relativistic case (Figure 3) the probability of R turns out to be almost independent from R. In other words, there will not be a majority of collisions with very small or very large energy exchange like in Figure 2, but each rate of energy redistribution has approximately the same probability.This marks a clear difference between low energy and high energy regimes.In Section 3 we shall use an agent-based simulation to see how this affects the equilibrium distribution in an analogous economic system.

Agent-Based Simulations and Saving Propensity
The analogies between the energy distribution in a gas and the income distribution in a society have been one of the main subjects of econophysics, since the pioneering works of the 1990s [9].In several works, agent-based simulations have been employed to investigate the relation between the rules of the microscopic interactions and the resulting equilibrium income distribution.A well established result is that a fixed saving propensity in the interactions gives rise to a distribution of the Gamma-kind, namely P(E) = kE β e −E/T , where E is here the income and T an effective temperature.This distribution does not have a flat tail.(P(E) is the probability density distribution of energy or money, in the econophysics analogy.The probability that an individual chosen at random has income in the infinitesimal interval (E, E + dE) is given by P(E)dE.The "temperature" T corresponds essentially, in the analogy, to the average wealth, whereby the Boltzmann constant is set by convention equal to 1.) The saving propensity, usually denoted with λ, is the fraction of its income that each individual "hedges" in the economic interactions, or in other words, individuals put at stake in the interactions a fraction (1 − λ) of their income.It is relatively straightforward to reproduce some of the mentioned result.For instance, Figure 4 shows the income distribution we have obtained from a simulation with 10 5 agents exchanging energy/money according to the rule [11] (i and j label two agents randomly chosen at each step; ξ is a random variable in the interval (0, 1)).
The system has been evolved to equilibrium with a series of 10 8 binary exchanges, and the whole procedure averaged for 50 times in order to minimize the effect of fluctuations.More in detail, we define, in a suitable programming language (Mathematica in our case) a vector E with 10 5 components.The component E i represents the wealth of the i-th agent of our idealized society.At the beginning all agents have the same wealth, e.g., E i = 20 units for any i in our simulation.Then the simulation algorithm choses at random two integer indices i, j between 1 and 10 5 .These are the agents who are currently interacting.Their wealths are changed according to Equation (15).The random choice of agents and the modification of their wealths is repeated for 10 8 times.This means that, on average, each agent has interacted with others for 10 3 times.(Thermalization times larger than 10 8 steps are found to give no changes in the final results.)A histogram of the wealth distribution is generated by rounding the wealth of each agent to the nearest integer and then using a counting function to define the frequency, for instance in Mathematica in our case with the code frequency = Table [Count[X,n],{n,0,120}], where X is the rounded income vector and 120 is a safe maximum value for the income.The frequency vector is the one whose components are shown in the histograms in Figures 4 and 5. Next we average the frequency vector over many realizations, in order to reduce the fluctuations, we normalize it by dividing each component by the sum of all components, and we use a FindFit function to fit it to the κ-distribution (16).This function employs a least-squares method and returns the best fit values of the three parameters in the κ-distribution.When the fit fails to show a power-law tail, it is because it returns negative values for α and β, which are not admissible for a κ-distribution.
A fit of the resulting distribution with the κ-distribution fails, and this confirms that no Pareto tail is present.Similar results can be obtained with other values of λ.Now let us run a simulation with microscopic exchange rules corresponding to the extreme relativistic case of Section 2.3.We set the energies of two particles after a collision as This means that, as found in Section 2.3, all energy exchange rates are equally probable, and as a consequence there is no fixed saving propensity.The result of a simulation with 10 5 agents, 10 8 thermalization steps and 100 averages is given in Figure 5.The curve can be fitted by a κ-distribution with a ratio α/κ 1.75, thus with a Pareto tail of exponent 2.75.Income distribution from an agent-based simulation (see details in text) with agents which exchange money in analogy with the behavior of relativistic particles, i.e., all exchange rates being equiprobable.The distribution has a Pareto tail proportional to 1/E 2.75 , where E is the income. of goods and services, while the money received may represent the profit from sales, a salary received etc.No distinction is made, at this level, between money as pure "cash" and other forms of wealth, like assets and properties, although in certain versions of the model such aspects have been considered, as well as the possibility that some individuals contract debts.This is clearly an over-simplification with respect to real economic interactions, but has the merit, from a complex systems perspective, to allow a micro-to-macro derivation of the income distribution starting from assumptions on the interactions between individuals, and not from the concept of rational representative individual and from the maximization of interest functionals like in classical economics theory.Boghosian [17,18] has extensively analyzed the mathematics of the Yard Sale model.In [19] he also gives an alternative economic interpretation: When two agents exchange a commodity, they may do so for a price that is different from the intrinsic value.This difference may be due to the agents' own personal utility functions-that is, their own senses of the value of the commodity.When this difference is nonzero, intrinsic wealth changes hands.

Income
There are further reasons to believe that kinetic models capture the "game-theoretical" strategies of the interacting agents much more than their physical origin could suggest.In fact, it was shown in [20] for game-theoretical models, that their large-numbers averages are equivalent to those of certain statistical mechanics systems.This suggests that some models of statistical mechanics originally developed to describe ensembles of inanimated atoms or spins are suitable also for the description of ensembles of individuals following a strategy, because in both cases the system tends to minimize some objective function.
A straightforward way of computing the average outcome of the elementary interactions in a wealth exchange model is to use an agent-based simulation [1,2].This is a tool widely employed in economics, also for models which comprise more sophisticated interactions.Each individual is defined as a variable in a programming code, which records its wealth and possibly other features.While executing a large number of discrete time-steps, an algorithm picks at each step two individuals at random and updates their wealth according to certain rules which simulate the interaction.A concrete example is given in Section 3 of this paper.
For a recent comparison of the results of an agent-based model with real data see [21].Saving propensity.In the simplest models of money transfer [9] the amount of money transferred at each encounter is typically independent from the money balances of the agents involved.The simulation usually proceeds by selecting two agents at random, choosing at random which of the two will be the payer and which the receiver, ad then transferring the fixed amount.This elementary step is repeated until the income distribution does not change.A different model, called the multiplicative asset exchange model, was introduced by Ispolatov et al. in 1998 [12].In this model the transferred amount is a fixed fraction γ of the payer's money.The stationary distribution of money in this model is similar, but not exactly equal, to the Gamma distribution density P(E) = cm β e −E/T .Using a Boltzmann kinetic equation, Ispolatov et al. have derived a relation between the parameters β and γ.Angle [22] connected the proportionality law heuristically introduced by Ispolatov et al. with the surplus theory of social stratification.Angle also considered a model where groups of agents have different values of γ, simulating the effect of education and other "human capital".Another model with an element of proportionality was proposed in Ref. [11].In this model the agents save some fraction λ of their money, whereas the rest becomes available for random exchanges.The exchange rules then become as in Equation ( 15) above.The coefficient λ is called the saving propensity and corresponds essentially to (1 − γ), where γ is the coefficient introduced in [12].

Figure 1 .
Figure 1.The function g(p 1x ) obtained by numerical solution of the differential Equation (12), for a strongly relativistic case (p/m = γv = 100 γ; v 1).We use units c = 1 throughout the paper, and in this plot we also take m = 1 for graphical reasons.

Figure 2 .Figure 3 .
Figure 2. Plot of the distribution of R in a Monte Carlo simulation of 5 × 10 6 almost classical collisions (p/m = γv = 0.01, γ 1, v 10 −2).R is expressed in % and is the fraction of kinetic energy transferred from Particle 1 to Particle 2. When p/m becomes even smaller, in the fully classical case, the peak for R → 0 (corresponding to a symmetric peak for R → 1) becomes higher.The dot at R = 51% is not meaningful.

Figure 4 .Figure 5 .
Figure 4. Income distribution from an agent-based simulation (see details in text) where agents exchange money with a fixed saving propensity.A fit with the κ-distribution confirms that there is no Pareto tail.