text:
# The pair approach applied to kinetics in restricted geometries: strengths and weaknesses of the method

## I Introduction

Classically, biochemical reaction kinetics is extrapolated from measurements in dilute solutions and fitted into the cellular reaction environment, and several flaws in this approach have been pointed out.  The main motivation of our work is to improve the understanding of diffusion-controlled reactions in topologically complex nanoscale environments represented in biological cells. In this study attempt is made to develop better theoretical methods which could describe diffusion-controlled reactions with boundaries. To achieve this goal, a particular way of doing calculation, the many-particle-density-function (MPDF) approach,  is modified to account for presence of boundaries.

The theoretical findings of this study are relevant for experimental work done in refs. . Even if we focus on biochemical reaction kinetics, the results should have an equal bearing on nanotechnological applications such as nanofluidics  or molecular electronics . Both are likely to be strongly dependent on reaction-diffusion behaviors of molecules (nanofluidics) or electrons and holes (molecular electronics) in restricted nanoscale geometries.

Most of the studies on diffusion-controlled reactions have been performed for infinite systems without boundaries and a variety of methods have been developed to do the analysis. The methods range from mean field treatments towards more exact approaches which employ quantum spin-chains , field theory , or MPDF formalism . The references  are an excellent review on the subject. The opposite case when reactions take place in restricted geometries with reactants confined into finite size, and eventually squeezed into very small volumes, is less understood. There is, however, some pioneering work in this area.  In here, focus is on testing performance of MPDF on diffusion-controlled reactions in finite volumes.

The infinite diffusion-controlled systems posses quite remarkable properties. When dimensionality of the system $`d`$ is lower than some critical dimension $`d_c`$ (e.g. for A+A reaction $`d_c=2`$ , and for $`A+B`$ $`d_c=4`$ ) a new non-trivial sort of kinetics sets in. Taking $`A+B0`$ as an example: Classical chemical kinetic rate equation for this reaction, with initial densities equal and homogeneous, is given by $`\dot{n}(t)=\lambda n(t)^2`$ where dot over symbol denotes time derivative. For large time $`t`$ this equation would predict density decay in the form of power law $`n(t)𝒜t^\alpha `$. The amplitude of decay equals $`𝒜=1/\lambda `$ and decay exponent is given by $`\alpha =1`$. In reality, $`\alpha =1`$ holds only for a sufficiently high dimensionality of the system when $`d>4`$, while for $`d<4`$ one has $`\alpha =d/4`$ and $`𝒜=\mathrm{const}\sqrt{n_0}D^{d/4}`$, and $`\mathrm{const}`$ is just a numerical factor. Please note that for $`d<4`$ decay amplitude $`𝒜`$ does not depend on the reaction rate $`\lambda `$ and exhibits lack of dependence on chemical details (universality).

The kinetics of the type described above is commonly referred to as anomalous or fluctuation-dominated. The term “anomalous” points to the fact that mean field (or classical rate equations) fail to describe such systems. The phrase “fluctuation-dominated” emphasizes importance of fluctuations in particle densities. Once the reaction creates a hole in the particle concentration, diffusion is very slow in restoring the homogeneous particle density. This has to do with recurrence of random walks. For $`d2`$ probability that the random walker will return to the same site after arbitrary number of steps is equal to one. Random walkers tend to wander around their initial position, and particles do not mix that well. Rule of thumb is that for lower dimensions the kinetics gets more anomalous. Role of dimensionality is well understood both for integer and fractal (non-integer) like dimensions.  On the other hand, much less is known what happens when one shrinks the system size, which is studied here.

To impact some progress in understanding reactions in restricted geometry we analyze performance of MPDF approach and modify it to account for finite reaction volume. To test such method of calculation the A+B model is used as a study case. The A+B model is natural choice for such a task. This model has been intensively studied for infinite system sizes.  It was found that the A+B reaction has the remarkable property that domains rich with A or B particles are formed as time goes on. Once domains are formed the reactions happen only at domain boundaries which leads to already mentioned decay exponent $`\alpha =d/4`$. It is interesting to study in what way the dynamics of the system (kinetics) changes as one reduces the volume available to reaction, in particular whether domain-like structure survives. For simplicity reasons, the one-dimensional (1d) case is studied. In the calculation that follows there is nothing special about 1d and present analysis can be easily extended to the two- or three- dimensional cases.

The A+B model in restricted geometry has been studied before with the assumption that one type of reactant is attached at the center of a small volume, and it was further assumed that one type of particle is in large excess.  The more realistic problem where all particles are allowed to move is much harder to solve, and the goal of our study is to describe such a situation. Also, in here, the focus is on the case when the initial number of reactants is the same, or roughly the same. Naturally, the shape of the reaction container might be important but this issue is not addressed at the moment. To avoid boundaries completely, our 1d system will be closed into the ring.

The paper is organized as follows. In section II the model is developed, i.e. detailed account is given of how particles move and react. Lattice model is used due to its conceptual simplicity. In section III equations of motion are derived using MPDF formalism in its simplest form (pair approach). In section IV equations of motion are solved analytically and it is shown how multi-exponential decay emerges. The results of computer experiments (Monte Carlo Simulations) are given in section V followed by a comparison between theory and Monte Carlo simulations in section VI. We conclude by analysis of strengths and weaknesses of the pair-approach applied to a reactions in restricted geometry in section VII. In appendix A details are given how to calculate effective reaction rate $`k(t)`$ which determines density decay. In appendix B the algorithm used to do Monte Carlo simulations is discussed in detail.

## II The lattice model

To test any theory one inevitably needs a model which serves as study case. The model used here is defined as follows. The two species, A and B, move on a 1d lattice performing random jumps with rates (diffusion constants) $`D_A`$ and $`D_B`$ respectively. It is assumed that $`D_A=D_B`$. Position of lattice sites is given by $`x_i=ih`$ with $`i=0,1,2,3,\mathrm{},M`$ and $`h`$ denotes lattice spacing. Sometimes, $`x`$ and $`y`$ will be used instead of $`x_i`$. Periodic boundary conditions are assumed and sites $`i=0`$ and $`i=M`$ are defined to be equivalent. There are M lattice sites in total and $`L=Mh`$. By using periodic boundary conditions it is possible to work with a system of finite size and yet keep the spatial translational invariance. This greatly facilitates the analytical treatment of the problem.

It is assumed that the reaction probability (per unit time) for particle $`A`$ at $`x`$ and $`B`$ at $`y`$ is given by $`\sigma (xy)`$. For $`\sigma (xy)`$ simplest possible form is used

$$\sigma (xy)=\sigma _0\theta (a|xy|).$$
(1)
where $`\theta (x)=0`$ for $`x<0`$ and $`\theta (x)=1`$ for $`x0`$. In this way two important aspects of chemical reactions are embedded, $`a`$ corresponds to the effective range of reaction and $`\sigma _0`$ is its strength. One could also say that each particle carries a ring of radius $`a/2`$ and when two rings overlap the particles can react. In this sense $`a/2`$ could be thought of as the size (radius) of particles. For simplicity reasons it is assumed that the reaction products do neither influence reactants, nor the A+B reaction. Also, exclusion or steric effects are not taken into account, i.e. particles are allowed to “enter” into each other (please see Fig. 1) and react with same probability independently from which direction they approach each other.

The model has the useful property that if $`a`$ is thought of as the size of reactants, then by varying $`a`$ several interesting situations can be studied. For example, when $`a`$ is on the order of the system size $`L`$, one can think of situations of extreme crowding. On the other hand when $`aL`$ reactants appear as point-like objects. In Fig. 1 we offer a schematic way how to think about these situations. The model presented above is solved analytically and numerically by a Monte Carlo simulation in the following sections.

## III Equations of motion in pair approximation

To solve the A+B reaction-diffusion model in a restricted geometry we use a many-particle-density-function formalism (MPDF), since it was already used to describe asymptotics of the same reaction in an infinite volume.  We modify the formalism and apply it to the case of a restricted geometry. In the following the formulation presented in ref  will be closely followed. On the way, the changes made to the original formalism will be discussed.

The dynamics of the system, as defined in previous section, is stochastic and governed by Master Equation which describes time evolution of configurational probabilities of the system $`P(c,t)`$,

$$\dot{P}(c,t)=\underset{c^{}}{}\left[W_{c^{}c}P(c^{},t)W_{cc^{}}P(c,t)\right]$$
(2)
where $`c`$ is short notation for occupancy of lattice sites and $`W_{cc^{}}`$ are transition probabilities which can easily be deduced from the previous description of the model. Here and throughout the paper dot over symbol denotes time derivative.

The quantities of interest are particle densities $`\rho _A(x,t)`$ and $`\rho _B(x,t)`$ and they can be calculated from $`P(c,t)`$ (at least in principle). Since system is closed into a ring translational invariance holds and concentrations cease to be position dependent which leads to $`\rho _A(x,t)=n_A(t)`$ and $`\rho _B(x,t)=n_B(t)`$. Following recipe in ref.  gives following equations for $`n_A`$ and $`n_B`$,

$`\dot{n}_A(t)=n_A(t)n_B(t){\displaystyle _{L/2}^{L/2}}𝑑x\sigma (x)Y(x,t)`$ (3)
$`\dot{n}_B(t)=n_A(t)n_B(t){\displaystyle _{L/2}^{L/2}}𝑑x\sigma (x)Y(x,t)`$ (4)
The $`Y(x,t)`$ denotes correlation function for AB pairs. Absence of correlations is signaled by $`Y(x,t)=1`$. Please note that in this work system size is finite which enters through finite integration domain in the integrals above (it might appear as minor technical detail but this fact is very important). Also, it is assumed that reversal symmetry holds, i.e. $`Y(x,t)=Y(x,t)`$.

Again, following ref.  one can derive equation for $`Y(x,t)`$ which is given by

$`\dot{Y}(x,t)=(D_A+D_B)Y^{\prime \prime }(x,t)\sigma (x)Y(x,t)`$ (5)
$`n_BY(x,t){\displaystyle _{L/2}^{L/2}}𝑑y\sigma (y)Y(y,t)\left[X_B(xy,t)1\right]`$ (6)
$`n_AY(x,t){\displaystyle _{L/2}^{L/2}}𝑑y\sigma (y)Y(y,t)\left[X_A(xy,t)1\right]`$ (7)
where prime denotes spatial derivative, and $`X_A(x,t)`$ and $`X_B(x,t)`$ correlation functions for AA and BB pairs respectively. The $`X_A(x,t)`$ and $`X_B(x,t)`$ obey similar equations which are not given here to save the space.

The equations (3)-(7) are derived under assumption of Kirkwood superposition approximation, which is a technical way of saying that dynamics is governed by pair effects. Naturally, assumption of the dominance of pairs effects is an approximation. It might or might not work, and the goal of present study is to test this. In the following, to make analytic treatment possible, equations will be simplified further by setting $`X_A(x,t)`$ and $`X_B(x,t)`$ equal to one. This amounts to ignoring correlations among AA and BB pairs. In ref.  it was shown that (for infinite reaction volume) such approximation is too severe and does not lead to correct decay exponent $`\alpha =d/4`$ (it gives $`\alpha =d/2`$). Nevertheless, in here we consider such simplification. The validity of such an approximation, together with the fact that we are using Kirkwood approximation, is tested via computer experiment later on.

The form of boundary conditions for $`Y(x,t)`$ differs from the one used in ref. . In the case of infinite system one takes

$$Y(x,t)1,x\mathrm{}$$
(8)
while for finite system with periodic boundary conditions another form has to be used

$$Y(x+L,t)=Y(x,t)$$
(9)
It will be shown later that the change from (8) to (9) leads to a qualitative change from power law to (multi) exponential behavior for correlation dynamics. The rest of the boundary conditions are standard, and are taken as in the case of an infinite system size,

$`n_A(0)=n_0`$ (10)
$`n_B(0)=n_0`$ (11)
$`Y(x,0)=1`$ (12)
Also, taking $`L\mathrm{}`$ should reproduce findings of ref. , within set of approximations employed here.

## IV Emergence of multi-exponential density decay

With assumptions $`X_A=X_B=1`$ Eq. (7) reduces to

$$\dot{Y}(x,t)=(D_A+D_B)Y^{\prime \prime }(x,t)\sigma (x)Y(x,t)$$
(13)
Eq. (13) is solved by using a Laplace transform as shown in the appendix A. To simplify the algebra, it is assumed that $`\sigma _0`$ is arbitrary large. The exponential behavior emerges due to the fact that the spectrum of Eq. (13) is discrete due to particular nature of the boundary conditions. The final expression for $`k(t)`$ reads

$$k(t)=k_{reg}(t)+2a\delta (t)$$
(14)
and details of calculation are given in the appendix A. The regular part of $`k(t)`$ is given by

$$k_{reg}(t)=\frac{8D}{L2a}\underset{m=1}{\overset{\mathrm{}}{}}e^{\kappa _mDt}$$
(15)
and $`\kappa _m`$ are constants of multi-exponential decay (eigenvalues),

$$\kappa _m=\pi ^2\left(\frac{2m1}{L2a}\right)^2$$
(16)
The $`\delta `$-function term in (14) arises from the second term on the right hand side of Eq. (A9) when $`\sigma _0\mathrm{}`$ (please see the appendix A).

Once $`k(t)`$ is available one can calculate $`n(t)`$ as

$$n(t)=\frac{n_0}{1+I(t)n_0+2an_0}$$
(17)
where $`I(t)=_0^t𝑑t^{}k_{reg}(t^{})`$ and

$$I(t)=\frac{8}{L2a}\underset{m=1}{\overset{\mathrm{}}{}}\frac{1}{\kappa _m}\left(1e^{\kappa _mDt}\right)$$
(18)
The $`2an_0`$ term in the denominator of (17) comes from the $`\delta (t)`$ term in Eq. (14). It describes the immediate annihilation of particles which are within reaction range. When $`\sigma _0\mathrm{}`$, this happens instantaneously. Thus there is a sudden jump in particle concentration. For finite $`\sigma _0`$ this jump becomes a smooth transition (exponential decay with decay exponent proportional to large number $`\sigma _0`$).

The question is whether one can obtain results for an infinite system from the Eq. (17) above. This can be done using Poisson resummation formula when $`Dt/(\frac{L}{2}a)^21`$. The Poisson resummation procedure gives $`I(t)t^{1/2}`$ which results in the wrong exponent for the density decay; $`n(t)t^{1/2}`$ instead of correct $`n(t)t^{1/4}`$. Thus we just reconfirm the well known fact that for infinite systems, the pair approach predicts too fast decay of particles. However, for finite systems, the situation is not that clear, it appears to depend on the type of initial conditions the real system is subject to.

In here we consider two types of initial conditions. (i) When initially there is an equal number of A and B particles, one has $`n0`$ as $`t\mathrm{}`$; and at the end all particles have to annihilate. (ii) One can look at an ensemble of similar systems with equal number of A and B particles at $`t=0`$ on average, $`N_{0,A}=N_{0,B}`$. In such a case, one has different asymptotics, $`N(t)N(\mathrm{})`$ as $`t\mathrm{}`$.

Theoretical prediction is that, as time goes to infinity, the particle density exponentially approaches the value $`n_{th}(\mathrm{})`$;

$$n_{th}(\mathrm{})=\frac{n_0}{1+n_0L}$$
(19)
The value for $`n_{th}(\mathrm{})`$ above can be obtained by sending $`t\mathrm{}`$ in Eqs. (17) and (18). From (19) one sees that asymptotically number of particles is given by

$$N_{th}(\mathrm{})=N_0/(1+N_0)$$
(20)
where $`N_{th}(\mathrm{})=Ln_{th}(\mathrm{})`$. Please note that $`N_{th}(\mathrm{})`$ never approaches zero and settles at a number between zero and one. In the case of type (i) initial conditions, all particles annihilate and $`N(\mathrm{})=0`$. This is clearly in contradiction with Eq. (19). However, situation is not that hopeless, as will be discussed later. For the type (ii) initial condition, for each member in ensemble there is a chance that more than one particle will be left, since one start dynamics with (random) excess of A or B particles at $`t=0`$. Thus, in average, $`N(\mathrm{})`$ will be larger than 0. Clearly, pair approach has more chance to describe this situation correctly.

In summary, we find an exponential decay in the long time limit which is a pure artifact of the finites of the system. There is a clear indication that the quality of prediction depends on the type of initial conditions used in experiment. Also, the approximations made in deriving Eq. (17) are rather severe and in order to check the applicability of such a pair-approach Monte Carlo simulation is used.

## V Results of Monte Carlo simulations of A+B reaction in restricted geometry

Figures 2, 3, and 4 summarize the results of the Monte Carlo simulations in $`d=1`$. The Monte Carlo algorithm is described in detail in appendix B. Figure 2 shows a simulation for a system with a large initial number of particles with a varying reaction range from a nearest-neighbor interaction with $`a/L=0.0001`$ towards a longer range with $`a/L=0.02`$. Figure 3 shows the case when there are initially very few (exactly 10=5A+5B) particles present in the reaction volume, also with varying reaction ranges from $`a/L=0.001`$ to $`a/L=0.2`$. Thus figures 2 and 3 give simulation results for type (i) initial condition. Figure 4 deals with type (ii) initial condition, when the initial number of particles in an ensemble varies with the constrain that the sum of A and B particles equals $`10`$. (For example, one run could be done with 7A and 3B, the other run with 5A and 5B, and a third run with 4A and 6B particles, etc.) Figure 1 is a sketch of how to think of various situations when $`a`$ changes from small to large values.

From figure 2 it can be seen that in the case of the nearest neighbor reaction range ($`a=1`$) four distinct regimes appear and the log-log plot is used to reveal them; (a) mean field decay, (b) the plateau region (c) power law decay and (d) exponential decay at the end. These regimes disappear as the reaction range is increased, and eventually, for very large $`a`$, one only has the exponential regime.

The mean field regime corresponds to annihilation of particles with all reactants being well mixed. This leads to depletion of lattice to concentration of the order $`n1/a`$, thus one particle per reaction range. Then diffusion starts to operate and mixes particles. What is interesting is that for very large values of $`a`$, the plateau region starts earlier and lasts longer. Apparently, it takes some time before the particles find each other by diffusion and start reacting again.

The power law decay starts after the plateau region. There is universality in the power-law regime since all curves with different values of $`a`$ merge into one. This is somewhat surprising since a larger $`a`$ should mean faster annihilation, which indeed happens in the mean field regime, but yet in the power law stage all curves share the same power-law behavior. We speculate that this has to do with self organization and build up of correlation.

The exponential regime is entered after the power law regime, when the number of particles in the system becomes small. With the present computer hardware it was not possible to resolve this exponential regime better. This is indeed done in figure 3 with a smaller lattice size and lower particle number.

To illuminate this exponential decay at later stages of annihilation, we performed simulations with a smaller number of particles (10=5A+5B) on a smaller lattice with $`10^3`$ sites. Thus we used type (i) initial condition. To obtain each curve we followed 1000-3000 realizations of dynamics and averaged over such an ensemble. The result is shown in figure 3. The upper figure is in log-$`t`$ scale to resolve the small and large $`t`$ region, respectively. The lower figure is in normal-$`t`$ scale and we use it to detect exponential decay (where a straight line indicates exponential decay).

The crossover from mean-field to plateau-like dynamics can be seen in the upper graph where all curves drop down to a plateau value which is $`a`$-dependent. The theoretical prediction for this plateau is $`n(0^+)=n_0/(1+2an_0)`$. The initial drop in concentration is large for large $`a`$-values. In the upper figure, it is hard to say when the plateau behavior turns into exponential decay.

The lower graph shows that decay is indeed exponential since density curves at late times are straight lines in the log$`{}_{10}{}^{}(n)`$-$`t`$ plot. Thus at the late times $`n\mathrm{exp}(\kappa _1t)`$. Also, the decay constant $`\kappa _1`$ is $`a`$-dependent since slopes are different for various values of $`a`$, and $`\kappa _1`$ becomes larger with larger $`a`$. Also, it appears that there is an upper limit for $`a`$ at which decay becomes infinitely fast. Naturally, this happens when $`a=L/2`$ since none of the particles can escape from each other. The qualitative dependence of $`\kappa _1`$ on $`a`$ just discussed is in agreement with theoretical prediction in Eq. (16) with $`m=1`$.

Figure 4 is obtained in a similar way as figure 2. The only difference is that figure 4 deals with the type (ii) initial condition. For the particular run, when $`N_{0,A}N_{0,B}`$, the final number of particles in the system is not zero. For example, when starting from 7A and 3B particles, the system will end up in the state of 4A particles. This comes from that fact that the A+B reaction conserves the particle difference $`N_A(t)N_B(t)=const`$. Curves for different values of $`a`$ saturate at one single value which is, independent of $`a`$. Clearly, the value of the plateau is solely controlled by the excess of particles at $`t=0`$, and can be calculated from theory if needed, but result of Monte Carlo simulation is equally informative.

## VI Comparison of computer experiment and theory

Figure 5 shows a comparison of the analytical treatment with computer experiment (simulation parameters as in figure 3). It can be seen that the pair (Smoluchowskii) approach does not predict that the number of particles in the system should approach zero. The reasons for this are discussed later but have to do with the fact that we are looking at highly symmetric situation with equal number of A and B particles all the time. To enforce such zero asymptotics by hand we use interpolation formula

$$n_{int}(t)=n(t)n(\mathrm{})(1e^{\kappa _1t})$$
(21)
where $`\kappa _1`$ is the first dominant large time exponent in expression for $`k(t)`$ (see Eqs. 17 and 18). It can be easily seen that the equation above holds exactly for $`t=0`$ and $`t=\mathrm{}`$. If Eq. (21) is used instead of (17) the agreement with simulation improves in the sense that the decay is exponential and the qualitatively theoretical exponent is roughly the same as the one obtained from simulations. More work is in progress to develop improved interpolation formulas.

Figure 6 deals with the same type of comparison, but with a simulation setup as in figure 4 when the initial number of particles is not fixed, just the total number (type (ii) initial condition). One can see that the agreement between theory and simulation is much better. Clearly, Smoluchowskii theory deals better with the type (ii) initial condition represented by the figure 4 than by the type (i) represented by the figure 3. Also, in figure 6, one can see that theory predicts too fast particle annihilation. This is no surprise since this is what one would expect form such pair approach which does not take into account formation of domains. (In that respect there is similarity with infinite systems, but only for the type (ii) initial condition).

## VII Discussion

The goal of present work was to impact some understanding of diffusion-controlled reactions in restricted geometries, with aim to describe some aspects of chemical reactions in biological cells. Two issues have been dealt with:

(1) The particular way of doing calculation was tested, the MPDF formalism developed by Kuzovkov and Kotomin (see ref.  for details). To be able to solve equations analytically the hierarchy of many-particle-densities was truncated at the level of three-particle-density using shortened Kirkwood superposition approximation. This approximation amounts to assuming that pair effects dominate correlations, and present calculation can be viewed as a variant of pair approach.

(2) A two species reaction-diffusion model $`A+B0`$ in a restricted geometry was taken as a study case. Two types of initial conditions were considered, type (i) where the initial number of A and B particles is strictly equal, and type (ii) initial conditions where the initial number of particles is equal only approximatively.

Thus the paper is best viewed as a method paper since the main goal is to test the strengths and weaknesses of the pair approach. To test quality of approximations involved all results have been compared with the results of computer experiment (Monte Carlo simulation).

From a theoretical point of view, it seems that the pair method, being widely used in calculation of bulk properties, works with mixed success for the restricted reaction diffusion systems, at least the one studied in here. The agreement between theoretical calculation and computer experiment is qualitative in the case of type (ii) initial conditions. In the case of type (i) initial conditions there is less agreement, however, situation is not that hopeless.

In the case of type (i) initial conditions pair approach makes error of the order of one particle (please see Eq. 20), since it predicts that, in the average, after very long time, there will be between zero and one particle in the system (though all particles should vanish). When initial number of particles is relatively large, the pair approach can describe evolution of system for rather long time, before the regime is reached where only one particle is left in the reaction volume. However, in the case of small initial number of particles, there is no such time interval, and the mismatch between pair approach and simulation has to be addressed more seriously.

The weakness of pair approach in dealing with type (i) initial condition rests on the fact that the truncated set of equations for many-particle-densities do not recognize any effects which go beyond pair correlations. For example, in the present case, all information related to the fact that initially there were 5A and 5B particles in the system, and that all particles have to vanish eventually, is missing. Work is in progress to pass such type of information from higher order particle-density-functions to lower order ones. For example, we already have better interpolation scheme than the one given in Eq. (21) for the case when there are initially three particles in the system, but we are trying to understand how to extend such analysis to higher numbers.

Interestingly enough, it seems that, contrary to the systems with infinite sizes, setting $`X_A(x,t)=X_B(x,t)=1`$ is reasonable approximation for a finite system, but we have to perform more tests. This could have to do with the fact that if the system is too small, there will be no time to develop clusters of A and B particles, and setting $`X_A(x,t)=1`$ and $`X_B(x,t)=1`$ might turn to be a good approximation after all. Thus, one does not have to turn to more complicated methods of calculation if qualitative results are needed. Nevertheless, it is highly desirable to see what happens as one includes correlations among AA and BB pairs.

To be able to solve equations analytically, we had to simplify MPDF considerably down to the level of a pair like approach and various calculation schemes contain pair approach as possible approximation. Perhaps the most common form of pair approach is the one suggested by Smoluchowskii (see e.g. ref.  for interesting review). The Smoluchowskii approach boils down to solution of Poisson equation with different boundary conditions. (Many readers will be familiar with this in the context of heat transfer or quantum mechanical problems.)

However, one has to keep in mind that pair approach is an approximation, and it has to be tested to see whether it works. (For example, the pair approach can not describe A+B reaction-diffusion model when system size is very large. This has been discussed in ref. .) The advantage of pair approach, in the form used here, is that it is possible to go beyond it in a systematic way.

Also, from the particular way we have approached pair problem, one can see that the difficulties associated with Eq. (20) are likely to be much deeper than just the fact that we are using pair approach. (This problem clearly vanishes when system size is infinite, as $`n_{th}(\mathrm{})`$ goes to zero.) Any scheme which focuses on low rank particle-density-function will suffer in a similar way, in the case of the highly symmetric, e.g. type (i), initial condition. One really has to find a way how to incorporate information of higher order correlation functions into lower order ones, without calculating higher order correlations functions explicitly. This is a pressing issue.

A few words about the model used. The goal of present work is to develop calculation method rather than to describe specific chemical system. We used the model which we could solve, within reasonable level of approximation. Nevertheless, the question whether present model has any relevance for real biological and chemical systems needs to be addressed.

The reaction diffusion model studied here appears to be too simple, for two reasons. First, it does not account for chemical details which enter only through two parameters, diffusion constant $`D`$ and reaction rate $`\lambda `$. For example, the exclusion effects and steric effects are not contained in it. Also, the influence of product molecule is completely ignored. Second, 1d character of the model might be too restrictive.

Despite its simplicity, the model used here contains basic characteristics of diffusion-controlled reactions (reaction times $`t_R`$ are much smaller than the corresponding diffusion times $`t_D`$); particles are moving on the on the lattice and react when within reaction range with no memory of initially velocity. Previous research reviewed in refs.  has shown that diffusion-controlled models, similar to the one used here, can be used to describe real chemical reactions. In particular, the A+B model have been used to study two reactions in capillary tube quite successfully, bromine + cyclohexene $``$ adduct, and Cu<sup>2+</sup> \+ disodium ethyl bis(5-tertrazolylazo)acetate trihydrate $``$ 1:1 complex in water.

Why can one be sloppy and ignore chemical details to some extent? The reason for this is universality. Most often, predictions of the reaction-diffusion models (on lattice) are insensitive to the details of the chemistry involved. (For example, the decay amplitude for A+B reaction does not depend on the reaction rate $`\lambda `$.) This statement is valid provided one deals with very large system sizes. Naturally, this view is not the only one. There are other ways to approach diffusion-controlled reactions. For more chemical or biological approach to diffusion-controlled reactions see refs.  and  respectively.

The simplicity of the model is not necessarily such a big handicap, until one reaches extremely small sizes. For small system sizes, density decay will start depending on details. However, there is a large window in system sizes between extremely large and extremely small where such kind of universality could survive. In here we push the model over its borders by studying situation of extreme crowding and not accounting for exclusion effects.

Furthermore, we would like to notify the reader that we do refer to the model here as a “toy” model, it relatively simple to formulate it. However, this does not mean that the model is easy to solve, quite the contrary. In the case of infinite system size it has taken a lot of research effort to clarify that the decay exponent indeed is $`d/4`$. This issue was finally settled in ref.  which provides a strict mathematical proof.

The model has a potential to describe experiments in refs.  where, for example, the average diameter of the reaction container (liposome) is $`L125\mu m`$. The reactants A (enzyme) and B (substrate) are of the size of $`a_E,a_S1nm`$ and the typical number of reactants inserted is on the order of $`N1000`$. Thus, $`aL`$ holds to a very good approximation. In these experiments reactants appear as point-like objects and there is no need to give structure to reactants. Also, problems associated with Eq. (20) will likely not to case any damage due to large number of particles at $`t=0`$.

From the data above, the concentration of particles $`n_{A,B}N/L^3`$, is easily estimated to be $`n1\mu m^3`$, and the typical distance between particles is $`d_{AB}1/c^{1/3}L/N^{1/3}0.1\mu m`$. Thus, in average, $`a_{S,E}d_{AB}`$, and particles are very well separated. Therefore, it is reasonable to assume that pair effects are the dominant ones. This in turn simplifies the theoretical description considerably. Clearly, there are other scales in the problem and the criteria on applicability of the pair approach are more subtle in reality.

To summarize, it would be extremely interesting to have a general method of calculation which could describe diffusion-controlled reactions in finite volumes, perhaps something on the level of the pair approach. Pair approach is attractive since inclusion of chemical details such as exclusion or steric effects is possible (see e.g. work in ref. ) which certainly opens a new route towards more quantitative results. However, pair approach is an approximation, and before burdening pair approach with increasing amount of chemical details, one has to test its ultimate reach. Present work is an attempt in this direction.

## A Derivation of the reaction rate $`k(t)`$

With initial condition $`Y(x,0)=1`$, the Laplace transform of (13) becomes

$$sY(x,s)1=(D_A+D_B)Y^{\prime \prime }(x,s)\sigma (x)Y(x,s)$$
(A1)
The correlation functions are symmetric around origin, i.e. $`Y(x,t)=Y(x,t)`$. They are also periodic in $`L`$. This implies that it is sufficient to focus on positive x axis and impose boundary conditions $`\frac{}{x}Y(x,s)=0`$ for $`x=0`$ and $`x=L/2`$. Equation (A1) is an ordinary second order differential equation, which is easily solved by solving it in regions $`0<x<a`$ and $`a<x<L/2`$ separately and then matching solutions at the end. After some algebra one obtains

$`Y_1(x,s)={\displaystyle \frac{1}{D}}({\displaystyle \frac{1}{\nu }}{\displaystyle \frac{1}{\mu }})\times `$ (A2)
$`\times {\displaystyle \frac{\mathrm{ch}(x\sqrt{\mu })}{\mathrm{ch}(a\sqrt{\mu })+\sqrt{\frac{\mu }{\nu }}\mathrm{sh}(a\sqrt{\mu })\mathrm{ch}[\sqrt{\nu }(\frac{\mathrm{L}}{2}\mathrm{a})]}}+{\displaystyle \frac{1}{D\mu }}`$ (A3)
and

$`Y_2(x,s)={\displaystyle \frac{1}{D}}({\displaystyle \frac{1}{\mu }}{\displaystyle \frac{1}{\nu }})\times `$ (A4)
$`\times {\displaystyle \frac{\mathrm{ch}[(\frac{L}{2}x)\sqrt{\nu }]}{\mathrm{ch}[(\frac{L}{2}a)\sqrt{\nu }]+\sqrt{\frac{\nu }{\mu }}\mathrm{sh}[(\frac{L}{2}a)\sqrt{\nu }]\mathrm{ch}(a\sqrt{\mu })}}+`$ (A5)
$`+{\displaystyle \frac{1}{D\nu }}`$ (A6)
where $`Y(x,s)=Y_1(x,s)`$ for $`0xa`$ and $`Y(x,s)=Y_2(x,s)`$ for $`axL/2`$ with $`\mu =(s+\sigma _0)/D`$ and $`\nu =s/D`$. The reaction rate $`k(s)`$ is given by $`k(s)=2\sigma _0_0^a𝑑xY_1(x,s)`$ and equals

$`k(s)`$ $`={\displaystyle \frac{2\sigma _0}{D}}({\displaystyle \frac{1}{\nu }}{\displaystyle \frac{1}{\mu }}){\displaystyle \frac{\mathrm{sh}(a\sqrt{\mu })}{\sqrt{\mu }}}\times `$ (A8)
$`\times {\displaystyle \frac{\mathrm{ch}(x\sqrt{\mu })}{\mathrm{ch}(a\sqrt{\mu })+\sqrt{\frac{\mu }{\nu }}\mathrm{sh}(a\sqrt{\mu })\mathrm{ch}[\sqrt{\nu }(\frac{L}{2}a)]}}+{\displaystyle \frac{2\sigma _0}{D\mu }}`$
We could not find the inverse Laplace transform of the expression above in closed analytic form. However, this is possible when $`\sigma _0\mathrm{}`$. In such a case one has

$$k(s)2\sqrt{\frac{D}{s}}\frac{1}{\mathrm{cth}[(L/2a)\sqrt{\frac{s}{D}}]}+\frac{2\sigma _0a}{\sigma _0+s}+𝒪(1/\sigma _0)$$
(A9)
The inverse Laplace transform of the approximate expression for $`k(s)`$ can be found by a residuum method. The $`s=0`$ is not a branching point nor pole. The only poles come from $`cosh`$ term in denominator which has poles at $`s_m=\pi ^2(m1/2)^2D/(L/2a)^2`$. This fully fixes form of $`k(t)`$ in Eq. (14).

## B Computer Experiment via Monte Carlo Simulations

We have chosen the minimal process algorithm for the simulations for two reasons. The first reason is that the algorithm reproduces the master equation (2).  Second, our goal is to study a whole range of particle sizes and relatively large numbers of particles at the same time. Clearly, there are another possibilities to carry out Monte Carlo simulation, but the main advantage of the minimal process algorithm is that it can be applied for systems containing relatively large number of particles. An original algorithm was devised for the situation where $`ah`$, i.e. particles react at the same lattice site or when nearest neighbors. We had to modify the original version of the algorithm to account for finite reaction range when $`ah`$. A detailed description of the algorithm is given bellow.

Algorithm:

* Site $`i`$ is chosen at random.
* If the site is empty go to step (5).
* For a chosen site $`i`$, one has to calculate the rate $`W_i`$ for a certain process to occur (diffusion or reaction). Also, one needs a null rate $`N_i`$ where nothing happens (the so called “null process”). The null rate is defined from $`W_i+N_i=Q`$, where $`Q`$ is arbitrary but known at each simulation step. $`Q`$ is chosen in such a way that none of the $`N_i`$ is negative. In practise, the case when $`Q`$ is taken as the largest of $`W_i`$ works best since this leads to the smallest possible values for $`N_i`$, i.e. chance that nothing is done in course of simulation is reduced. (Please note that this requires that $`Q`$ is updated as simulation proceeds, but can be done in a straight forward manner as explained in ref. )

$`W_i=D_i+R_i`$ accounts for possibilities that a particle at the site diffuses to the neighboring site with rate $`D_i`$, or reacts with a particle in some other site with rate $`R_i=_{j\mathrm{\Omega }_i}\sigma (r_{ij})`$. $`\mathrm{\Omega }_i`$ denotes set of sites which are within reaction range of the site $`i`$. The calculation of $`R_i`$ is by far the most costly step when $`a`$ is large. In that case, a large region has to be searched in order to find all particles within $`\mathrm{\Omega }_i`$. This step costs $`M_{search}(a/h)^d`$ computational steps if the sites are checked one by one. The cost can be reduced further by introducing a list which specifies which sites that contain particles within the reaction range of the particle at site $`i`$. In that case one has to update the list for each diffusion step made. The best algorithm we have so far updates the list in roughly $`M_{serach}(a/h)^{d1}`$ steps.
* Once the rates for the site $`i`$ have been calculated one can use them to evaluate probabilities for specific process $`p^{(D)}=D_i/Q`$, $`p^{(R)}=R_i/Q`$ and $`p^{(null)}=N/Q`$. Once the probabilities are calculated a certain process is chosen by linear selection algorithm. First one decides if diffusion, reaction or nothing is going to happen. If diffusion is to happen then the particle is moved to one of the randomly chosen $`2\times d`$ nearest neighbors. If reaction was chosen, then one of the sites containing particles in reaction range is chosen at random, e.g. at site $`j`$, and pair of particles from site $`i`$ and $`j`$ are annihilated.
* Time is updated according to the formula $`tt+\mathrm{\Delta }t`$ where $`\mathrm{\Delta }t=1/LQ`$ where $`L`$ was specified before and $`Q`$ is the maximum rate at the present step.
* Move back to (1) unless some criteria to stop is invoked.

Applying the same type of reasoning as in ref.  one can see that the algorithm proposed here reproduces the behavior described by the master equation (2). As time of the simulation progresses we monitor the number of particles and calculate all the statistics.

As the original minimal process algorithm, the present simulation method is not that efficient at the later stages of dynamics when the lattice becomes sparse. The quantity that governs computational cost of this method is the number of Monte Carlo steps needed to see some change in the number of particles. We describe it by the number of Monte Carlo steps needed to annihilate the last pair of particles.

To make such estimate, it is best to move to the reference frame of one of these particles. Then one particle is fixed and another one is trying to find it. The number of diffusion steps that the moving particle needs to find the one who sits still is roughly given by $`M(\mathrm{diff})L^d/a`$ (here and in the following it is implicitly assumed that every length variable is measured in units of lattice spacing $`h`$). Each diffusion step bears $`M(\mathrm{step}/\mathrm{diff})`$ computational steps which gives the total number of steps to annihilate the pair of particles equal to $`M(\mathrm{tot})=M(\mathrm{diff})M(\mathrm{step}/\mathrm{diff})`$. The number of Monte Carlo steps per one diffusion step is roughly 1, $`M(\mathrm{step}/\mathrm{diff})1`$. However, calculation of $`R_i`$ requires updating the internal list which costs $`M_{search}a^{d1}`$ search steps whenever the particle is moved. Thus, the true number of computational steps per diffusion step is given by $`M(\mathrm{step}/\mathrm{diff})M_{search}a^{d1}`$. Finally, one gets an estimate for the number of computational steps needed to annihilate the last pair of particles as $`M(\mathrm{tot})L^da^{d2}`$.

The algorithm has an interesting property that for $`d=1`$ there is a reduction in the computational cost when comparing large and small $`a`$ cases. For larger $`a`$ the algorithm works more efficiently. For $`d=2`$ the computational cost does not depend on $`a`$. Simulating a large $`a`$ situation for $`d=3`$ is more costly. One could avoid this growing cost problem at $`d=3`$ by browsing through particles instead of searching for sites when calculating $`R_i`$. This is clearly the preferred option when the number of particles in the system is not that large.


format: text

text:
# Saturation of optimality limits in hadron-hadron scatterings.

## Abstract

In this paper the optimal unitarity lower bound on logarithmic slope of the diffraction peak is investigated. It is shown that the unitarity lower bound is just the optimal logarithmic slopes predicted by the principle of least distance in space of states. A systematic tendency towards the saturation of the most forward-peaked limit is observed from all the available experimental data of all principal hadron-hadron (e.g., $`PP,`$ $`\overline{P}P,`$ $`K^\pm P,`$ $`\pi ^\pm P`$ ) scatterings practically at all laboratory momenta.

1. Introduction

Recently, in Ref. , by using reproducing kernel Hilbert space (RKHS} methods \[2-4\], we described the quantum scattering of the spinless particles by a principle of minimum distance  in the space of the scattering states (PMD-SS). Some preliminary experimental tests of the PMD-SS, even in the crude form  when the complications due to the particle spins are neglected, showed that the actual experimental data for the differential cross sections of all $`PP,`$ $`\overline{P}P,`$ $`K^\pm P,`$ $`\pi ^\pm P,`$ scatterings at all energies higher than 2 GeV, can be well systematized by PMD-SS predictions. Moreover, connections between the optimal states , the PMD-SS in the space of quantum states and the maximum entropy principle for the statistics of the scattering channels was also recently established by introducing quantum scattering entropies \[5-7\]. Then, it was shown that the experimental pion-nucleon as well as pion-nucleus scattering entropies are well described by optimal entropies and that the experimental data are consistent with the principle of minimum distance in the space of scattering states (PMD-SS) .

An important model independent result obtained Ref. , via the description of quantum scattering by the principle of minimum distance in space of states (PMD-SS), is the following optimal lower bound on logarithmic slope of the forward diffraction peak in hadron-hadron elastic scattering:

$$b\frac{d}{dt}\mathrm{ln}\left[\frac{d\sigma }{d\mathrm{\Omega }}(s,t)\right]_{t=0}b_o\frac{\overline{\lambda }^2}{4}[\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(1)1]$$
(1)
where$`\sqrt{s}`$ and $`\sqrt{\left|t\right|}`$ are c.m. energy and transfer momentum, respectively, while$`\overline{\lambda }`$ is c.m. d’ Broglie wave length, $`\frac{d\sigma }{d\mathrm{\Omega }}(x)`$ is the differential cross section, $`x\mathrm{cos}\theta ,`$ being the c.m. scattering angle, and $`\sigma _{el}`$ is the integrated elastic cross section. So, the quantum scattering state which saturate the optimal bound (1) is the most forward-peaked quantum state. The optimal bound (1) improves in a more general and exact form not only the unitarity bounds derived by Martin , MacDowell and Martin  (see also Ref. ) for the logarithmic slope $`b_A`$ of absorptive contribution $`\frac{d\sigma _A}{d\mathrm{\Omega }}(s,t)`$ to the elastic differential cross sections but also the unitarity lower bound derived in Ref.  for the slope $`b`$ of the entire $`\frac{d\sigma }{d\mathrm{\Omega }}(s,t)`$ differential cross section. Therefore, it would be very important to clarify the optimal status of the bound (1) and to make an experimental detailed investigation of its saturation in the hadron-hadron scattering especially in the low energy region.

The aim of this paper is twofold: to show that the optimal bound (1) is the ”singular” solution of the problem to find an unitarity lower bound on the logarithmic slope b when: $`\sigma _{el}`$ and $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ are given, and, to obtain experimental tests of the bound (1) at all available energies for the principal hadron-hadron elastic scatterings: $`\pi ^\pm P\pi ^\pm P,`$ $`K^\pm PK^\pm P,`$ $`PPPP,`$ and $`\overline{P}P\overline{P}P`$. So, in Sect. 2 the problem of the optimal status of the bound (1) is completely solved by showing that this bound is the singular solution ($`\lambda _0=0`$) in the unconstrained Lagrangean minimization problem. The relation of the optimal bound (1) with the Martin-MacDowell bound as well as its saturations in different particular cases (e. g. pure absorptive cases, etc.) are discussed. Some details about the PMD-SS-history as well as about other predictions of the principle of minimum distance in the space of states are presented in Sect. 3. Then, it is shown that the optimal state obtained with the simplified version of the PMD-SS have not only the property that is the most forward-peaked quantum state but also possesses many other peculiar properties such as: (i) scaling in the variable $`\tau _o=2\sqrt{|t|b_o}`$ , (ii) maximum Tsallis-like scattering entropies  $`S_\theta (q),`$ $`S_L(q)`$, $`S_{\theta L}(q)`$, $`\stackrel{~}{S}_{\theta L}(q)`$ for the all positive nonextensivity parameters $`q`$, etc., that make it a good candidate for the description of the quantum scattering. The results on the experimental test of the saturation of the most forward-peaked state limit are presented in Sect. 4 while the Sect. 5 is reserved for summary and conclusions.

2. The problem of the unitarity lower bound on logarithmic slope

Let us start with the hadron-hadron elastic scattering process of spinless particles for which the description of the scattering amplitude $`f(x)`$ of the system is given in terms of partial amplitudes $`f_l,l=0,1,\mathrm{}`$as

$$f(x)=(2l+1)f_lP_l(x),x[1,1],f_lC$$
(2)
where $`P_l(x),l=0,1,\mathrm{}`$ , are Legendre polynomials, $`f_l,l=0,1,\mathrm{}`$ being the partial amplitudes. If the normalization is chosen such that $`\frac{d\sigma }{d\mathrm{\Omega }}(x)=f(x)^2,`$ the elastic integrated cross section $`\sigma _{el}`$, is expressed in terms of partial amplitudes by

$$\frac{\sigma _{el}}{2\pi }=2(2l+1)f_l^2=f^2$$
(3)
Now, a rigorous unitarity lower bound on the slope parameter b can be obtained by solving completely the following constrained minimization problem:

$$\mathrm{min}\{b\},when\sigma _{el}=fixed,\frac{d\sigma }{d\mathrm{\Omega }}(1)=fixed,f_l^2Imf_l,l=0,1,..$$
(4)
which is equivalently to the following unconstrained minimization problem:

$$L\lambda _0b+\lambda _1\left[\frac{\sigma _{el}}{4\pi }(2l+1)f_l^2\right]+\lambda _2\left[\frac{d\sigma }{d\mathrm{\Omega }}(1)(2l+1)f_l^2\right]+\xi _lg_l\mathrm{min}$$
(5)
where $`g_lImf_lf_l^2,`$ $`\lambda _i,i=0,1,2,`$ and $`\xi _l`$ , $`l=0,1,\mathrm{},`$ are the corresponding Lagrange multipliers, and $`b=(R\stackrel{}{R}+I\stackrel{}{I})/(R^2+I^2),`$ where $`R(2l+1)Ref_l,`$ $`I(2l+1)Imf_l`$, and $`\stackrel{}{R}(2l+1)\frac{l(l+1)}{2}Ref_l`$, $`\stackrel{}{I}(2l+1)\frac{l(l+1)}{2}Imf_l`$.

The non singular solution $`(\lambda _00)`$ of the minimization problem (5) was completely treated in Ref. . Here we must underline that, according to the rules from Ref. , the singular solution $`\lambda _0=0,`$ if it exists, can be essential in obtaining a rigorous unitarity lower bound on the logarithmic slope parameter $`b`$. In order to prove this statement, let us consider as an simple example the minimization problem: $`x\mathrm{min}`$, subject to $`x^2+y^2=0`$. The solution of this problem is evident: $`x_{\mathrm{min}}=y_{\mathrm{min}}=0`$. Now, let use the Lagrange multiplier’s problem with $`\lambda _00`$: $`L=x+\lambda (x^2+y^2).`$ Then, we obtain the following incompatible equation: $`1+2\lambda x=0`$, $`2\lambda y=0`$, and $`x^2+y^2=0`$. Hence, we must consider the case $`\lambda _0=0`$ from which we get the correct solution. Such a singular solution of the Lagrange function (5) can be obtained by solving the minimization problem

$$\left\{(2l+1)f_l^2+\alpha \left[\frac{d\sigma }{d\mathrm{\Omega }}(1)(2l+1)f_lP_l(1)^2\right]\right\}\mathrm{min}$$
(6)
The problem (6) is just the problem of minimum norm (minimum distance in the space of states) when $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ is fixed. Thus, the unique solution of (6) exists and is given the following optimal state (see Ref. )

$$f_{o1}(x)=f(1)\frac{K(x,1)}{K(1,1)}=f(1)\frac{\underset{L_o+1}{\overset{.}{P}}(x)+\underset{L_o}{\overset{.}{P}}(x)}{\left(L_o+1\right)^2}$$
(7)
where $`\underset{l}{\overset{.}{P}}(x)dP_l(x)/dx`$while the reproducing kernel $`K(x,1)`$ is given as follows

$`K(x,1)=\frac{1}{2}_{l=0}^{L_o}(2l+1)P_l(x)=\frac{1}{2}\{\underset{L_o+1}{\overset{.}{P}}(x)+\underset{L_o}{\overset{.}{P}}(x)\}`$,$`K(1,1)=\left(L_o+1\right)^2/2`$,

and

$$L_o=integer\left\{\left[\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(1)\right]^{1/2}1\right\}$$
(8)
Hence, the singular solution ($`\lambda _0=0,\chi _l=0,l=0,1,\mathrm{}.`$) of the minimization problem (5) (or (4)) is just the optimal slope

$$b_o=\frac{d}{dt}\mathrm{ln}[f_{o1}(x)^2]_{t=0}=\overline{\lambda }^2\frac{L_o(L_o+2)}{4}=\frac{\overline{\lambda }^2}{4}[\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(1)1]$$
(9)
This result is in agreement with that obtained in Ref.  directly from the fundamental inequality $`f(y)^2\frac{\sigma _{el}}{2\pi }K(y,y),`$ $`y[1,+1].`$ Then, by developing the scattering amplitude $`f(y)`$ and the reproducing kernel function $`K(y,y)`$ in Taylor’s series around the point $`y=1`$ and $`L=L_o`$, it was shown that

$$\frac{d}{dy}\left[\frac{d\sigma }{d\mathrm{\Omega }}(y)\right]_{y=1}\frac{\sigma _{el}}{2\pi }\frac{d}{dy}\left[K(y,y)\right]_{y=1,L=L_o}=\frac{d\sigma }{d\mathrm{\Omega }}(1)\frac{L_o(L_o+2)}{2}$$
(10)
From (10) and (8) we get the optimal lower bound (1) since, according to the definition of the logarithmic slope $`b`$ combined with the kinematical relations $`t=2q^2(1y)`$, we have $`dt=2dy/\overline{\lambda }^2`$ and

$$b=\frac{\overline{\lambda }^2}{2}\frac{d}{dy}\left[\frac{d\sigma }{d\mathrm{\Omega }}(y)\right]_{y=1}/\frac{d\sigma }{d\mathrm{\Omega }}(1)$$
(11)
Now, by comparing the singular solution (1) with that $`b_{\mathrm{min}}(\lambda _00)`$ obtained in Ref.  for the non singular case $`\lambda _00,`$ we see that: $`b_o>b_{\mathrm{min}}(\lambda _00)\frac{8}{9}b_o`$. Consequently, since we get the chain of bounds: $`bb_o>b_{\mathrm{min}}(\lambda _00),`$ the unique solution of the problem (5) is that given by the optimal bound (1) since max$`\{b_{\mathrm{min}}(\lambda _0=0),b_{\mathrm{min}}(\lambda _00)\}=b_o`$

Remark 1: The solution to the minimization problem (5) is of relevance since in the pure absorptive case ($`Ref(x)=0`$) the singular solution $`bb_o`$ as well as that non singular $`b>b_{\mathrm{min}}(\lambda _00)\frac{8}{9}b_o`$ (see Ref.) are going in

$$b_A\frac{\overline{\lambda }^2}{4}\left[\frac{\sigma _T^2}{4\pi \overline{\lambda }^2\sigma _{el}}1\right]b_o^A$$
(12)
and, respectively, in the MacDowell-Martin bound

$$b_A>\frac{8}{9}b_o^A=\frac{2\overline{\lambda }^2}{9}\left[\frac{\sigma _T^2}{4\pi \overline{\lambda }^2\sigma _{el}}1\right]b_{MDM}$$
(13)
since

$$\frac{d\sigma ^A}{d\mathrm{\Omega }}(1)=[Imf(x)]^2=\sigma _T^2/(4\pi \overline{\lambda })^2,\mathrm{andd}b=b_A$$
(14)
Remark 2: Let $`\sigma _{el}`$ and $`\sigma _T`$ be fixed from experiment let us consider the minimization problem

$$\mathrm{min}\{b\},\mathrm{when}\sigma _{el}=fixed,\mathrm{and}\sigma _T=fixed,f_l^2Imf_l,l=0,1,.$$
(15)
which is equivalent to the minimization of the Lagrangean

$$L\lambda _0b+\lambda _1\left[\frac{\sigma _{el}}{4\pi }(2l+1)f_l^2\right]+\lambda _2\left[\sigma _T(2l+1)Imf_l\right]+\xi _lg_l\mathrm{min}$$
(16)
Then, the unique solution of Eq. (15) is given by

$$b\mathrm{max}\{b_{\mathrm{min}}(\lambda _0=0),b_{\mathrm{min}}(\lambda _00)\}=\mathrm{max}\{b_o^A,b_{MDM}\}=b_o^A$$
(17)
The quantum scattering state which saturate the optimal bound (23) is the most forward-peaked quantum state which is obtained as the absorptive limit of the optimal state (7).

Indeed, the singular solution $`\lambda _0=0`$ of the problem (16), can be obtained by solving the minimization problem

$$\left\{(2l+1)f_l^2+\alpha \left[\sigma _T(2l+1)Imf_l\right]\right\}\mathrm{min}$$
(18)
The unique solution of this problem is given by

$$\begin{array}{c}Ref_l^o=0,\mathrm{forall}l,\\ Imf_l^o=4\pi \overline{\lambda }\frac{\sigma _{el}}{\sigma _T},\mathrm{for}\mathrm{\hspace{0.33em}0}lL_o^A,Imf_l=0,\mathrm{for}lL_o^A+1\end{array}$$
(19)
where

$$L_o^A=\mathrm{integer}\left\{\left[\frac{\sigma _T^2}{4\pi \overline{\lambda }^2\sigma _{el}}\right]^{1/2}1\right\}=\mathrm{integer}\left\{\left[\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}^A(1)\right]^{1/2}1\right\}$$
(20)
By introducing the results (19) and (20) in Eq. (2) we obtain

$$\begin{array}{c}Ref_{o1}^A(x)=0\\ Imf_{o1}^A(x)=Imf(1)\frac{{\scriptscriptstyle (2l+1)P_l(x)}}{(L_o^A+1)^2}=Imf(1)\frac{\underset{L_o+1}{\overset{.}{P}}(x)+\underset{L_{oT}}{\overset{.}{P}}(x)}{\left(L_o^A+1\right)^2}\end{array}$$
(21)
Hence, for purely absorptive amplitudes the optimal solution to the minimization problem (19) exists and is given by Eq. (21) which in fact is the pure absorptive limit of the optimal state (7) but with $`L_o=L_o^A.`$ Hence, in the pure absorptive case from (20) we get

$$b_{\mathrm{min}}(\lambda =0)=\overline{\lambda }^2\frac{L_o^A(L_o^A+2)}{4}=\frac{\overline{\lambda }^2}{4}\left[\frac{\sigma _T^2}{4\pi \overline{\lambda }^2\sigma _{el}}1\right]b_o^A$$
(22)
Therefore, the unique solution of the problem (15) is given by

$$b\mathrm{max}\{b_o^A,b_{MDM}^A\}=b_o^A$$
(23)
The MacDowell-Martin bound $`b^A>b_{MDM}^A`$ (13) emerge as the non singular solution in a special case of the minimization problem (5) when $`Ref(x)=0`$ and due the fact that is $`\frac{8}{9}b_o^A`$ is eliminated in the final solution (23) of the problem (15). Hence, the bound (1) include in a more general and exact form the MacDowell-Martin bound since the validity of the bound $`bb_o`$ includes the validity of any inequality obtained as a consequence of the inequality $`b_o>a`$, e.g., $`ab_{MDM}^A`$(13).

Finally, we note that the fundamental inequality (16) from Ref.  tell us that if $`\sigma _{el}`$ and$`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ (or if $`\sigma _{el}`$ and $`\sigma _T`$) are fixed from experiment, then the number $`(L+1)`$ of partial amplitudes in any phase shift analysis (PSA) must obey the optimal bounds

$$\left(L+1\right)^2\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(1),\left(L+1\right)^2\frac{\sigma _T^2}{4\pi \overline{\lambda }^2\sigma _{el}}$$
(24)
respectively, or equivalently $`LL_oL_o^A.`$

3. Why principle of minimum distance in space of states?

As we already underlined in introduction, in Ref.  are investigated the essential features of the hadron-hadron scattering by using the optimization theory  in the Hilbert space of the partial scattering amplitudes. Then, knowledge about the hadron-hadron scattering system (or more concretely, about partial amplitudes) are deduced by assuming that it behaves as to optimize some given measure of its effectiveness, and thus the behavior of the system is completely specified by identifying the criterion of effectiveness and applying constrained optimization to it. This approach is in fact known as describing the system in terms of an optimum principle.

The earliest optimum principle was proposed by Hero of Alexandria (125 B.C.) in his Catoptrics in connection with the behavior of light. Thus, Hero of Alexandria (125 B.C.) proved mathematically the following first genuine scientific minimum principle of physics.

Hero’s principle of minimum distance (PMD) tell us that when a ray of light is reflected by a mirror, the path actually taken from the object to the observer‘s eye is the shortest path from all possible paths. If we extend the PMD-idea to the behavior of light in gravitational fields, then, we obtain immediately that according to the PMD-optimum principle, modified to include the interaction of light with the gravitational field, the light must move on a specific shortest path which is the geodesic. Having in mind this successful result, the principle of minimum distance was recently  extended to the quantum physics by choosing the ”partial transition amplitudes” as fundamental physical quantities since they are labelled with good quantum numbers such as charge, angular momentum, isospin, etc. These physical quantities are chosen as the system variational variables  while the distance in the Hilbert space of the quantum states is taken as measure of the system effectiveness expressed in terms of the system variables. The principle of minimum distance in the space of states (PMD-SS) is chosen as variational optimum principle by which one should obtain those values of the ”partial amplitudes” yielding minimum effectiveness. Of course this new optimum principle can be formulated in a more general mathematical form by using the S-matrix theory of the strong interacting systems.

The Hilbert space H of the scattering states, with the inner product $`<.`$$`,.>`$ and the norm $``$, defined on the interval $`S(1,+1)`$, can be given by

$$<f,g>=_1^{+1}f(x)\overline{g(x)}𝑑x,f,g𝐇,$$
(25)
$$f^2=<f,f>=_1^{+1}f(x)^2dx,f𝐇$$
(26)
If the quantum distance $`D(f,g)`$ between two quantum scattering states $`f`$ and $`g`$ is defined  by

$$D(f,g)=\underset{\mathrm{\Phi }}{\mathrm{min}}fg\mathrm{exp}(i\mathrm{\Phi })=[f^2+g^22<f,g>]^{\frac{1}{2}}$$
(27)
then the PMD-SS can be formulated in the following form.

PMD-SS: If $`h`$ is the quantum state of the system when the interaction is missing, then the true interacting quantum state $`f`$ of the interacting system is that state which have the shortest distance $`D(f,h)`$ in the space of interacting states compatible with the constraints imposed by the interaction. (Of course, if $`h=0`$, then $`D(f,0)=f`$.)

Therefore, let $`M_{eff}`$ be the measure of the scattering effectiveness and let us consider the problems of optimization of few $`M_{eff}`$ defined as follows:

(a) $`M_{eff}D^2(f,0),`$ or equivalently $`M_{eff}=\sigma _{el}/4\pi `$, where $`\sigma _{el}`$ is given in terms of partial amplitudes by Eq. (3);

(b) $`M_{eff}b`$, where the logarithmic slope $`b`$ of the diffraction peak is expressed in terms of the partial amplitudes as in

Sect. II;

(c) $`M_{eff}S_L=(2l+1)p_l\mathrm{ln}p_l,`$ where: $`p_l=4\pi |f_l|^2/\sigma _{el},`$ and $`(2l+1)p_l=1;`$

(d)$`M_{eff}S_\theta `$ $`=_1^{+1}𝑑xP(x)\mathrm{ln}P(x)`$, where $`P(x)=\frac{2\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(x)`$, with $`_1^{+1}𝑑xP(x)=1;`$

Now, let us consider the constrained optimization problems:

$$\mathrm{extremum}\{M_{eff}\}\mathrm{when}\sigma _{el}\mathrm{and}\frac{d\sigma }{d\mathrm{\Omega }}(1)\mathrm{are}\mathrm{fixed}.$$
(28)
which, via Lagrange multipliers method , are reduced to the problems:

$$\mathrm{\pounds }_F\lambda _0M_{eff}+\lambda _1\left[\frac{\sigma _{el}}{4\pi }(2l+1)f_l^2\right]+\lambda _2\left[\frac{d\sigma }{d\mathrm{\Omega }}(1)(2l+1)f_l^2\right]\mathrm{extremum}$$
(29)
Hence, according to the general Lagrange multiplier rules , for all these optimization problems (28)-(29) we are obligated to analyze in detail the singular solution $`\lambda _0=0`$. Such a singular solution of the Lagrangean function (29) can be obtained by solving the minimization problem (6) which is just the problem of minimum norm (minimum distance in the space of states)

$$\mathrm{min}[D(f,0)]\mathrm{when}\frac{d\sigma }{d\mathrm{\Omega }}(1)\mathrm{is}\mathrm{fixed}fixed$$
(30)
Thus, the unique solution of the problem (30) exists and is given by the PMD-SS-optimal state (7) (see details in Ref. ). Therefore, the optimal effectiveness are as follows (see e.g. Ref.  for the (c), (d) cases):

$$\begin{array}{c}M_{eff}^{o1}=b_o,\mathrm{in}\mathrm{case}(\mathrm{b})\\ M_{eff}^{o1}=S_L^{o1}=\mathrm{ln}[2K(1,1)=\mathrm{ln}[(L_o+1)^2],\mathrm{in}\mathrm{case}(\mathrm{c})\\ M_{eff}^{o1}=S_\theta ^{o1}=_1^{+1}𝑑x\left[\frac{[K(x,1)]^2}{K(1,1)}\right]\mathrm{ln}\left[\frac{[K(x,1)]^2}{K(1,1)}\right],\mathrm{in}\mathrm{case}(\mathrm{d})\end{array}$$
(31)
while the optimal angular momentum is given by Eq. (8).

Now, we recall briefly some peculiar properties of the PMD-SS-optimal state (7) obtained in Refs. , - that can illustrate the importance of the PMD-SS for an unified description of all quantum phenomena. These characteristic features are as follows:

I. The PMD-SS-optimal scaling property (see the experimental tests in Fig. 1 in Ref.)

$$\frac{1}{\frac{d\sigma }{d\mathrm{\Omega }}(s,1)}\frac{d\sigma ^{o1}}{d\mathrm{\Omega }}(s,x)=\left[\frac{K(x,1)}{K(1,1)}\right]^2\left[\frac{2J_1(\tau _o)}{\tau _o}\right]^2,\mathrm{for}\mathrm{small}\tau _o$$
(32)
as a function of the PMD-SS-optimal scaling variable $`\tau _o`$

$$\tau _o2\left[\left|t\right|b_{o1}\right]^{1/2}=\left\{\overline{\lambda }^2\right|t|\left[\frac{4\pi }{\sigma _{el}}\frac{d\sigma }{d\mathrm{\Omega }}(1)1\right]\}^{1/2}$$
(33)
II. The PMD-SS-optimal state (7) is the most forward peaked scattering state in the sense given by inequality (1) (see Figs. 1-3 in this paper);

III. The PMD-SS-optimal state (7) is the state with maximum scattering entropies (see proofs and the experimental illustrations in Ref. ) in the sense that the entropies $`S_L`$and $`S_\theta `$ of any scattering state must fulfill the inequalities

$$S_\theta S_\theta ^{o1},S_LS_L^{o1}$$
(34)
IV. PMD-SS-optimal state (7) is the state of equilibrium of the angular momenta channels considered as a quantum statistical ensemble.  Indeed, using Eqs. (31) and (34) we observe that the entropy $`S_L`$ is similar to the Boltzmann entropy with a maximum value given by the logarithm of number of the optimal states: $`S_LS_L^{o1}=\mathrm{ln}[(L_o+1)^2].`$The results like (34) are also proved in Ref.  for the nonextensive statistics  of the angular momenta channels.

V. From mathematical point of view, the PMD-SS-optimal states, are functions of minimum constrained norm and consequently can be completely described by reproducing kernel functions (see Ref. \[1,3-4\].  So, with this respect the PMD-SS-optimal states from the reproducing kernel Hilbert space (RKHS) of the scattering amplitudes are analogous to the coherent states from the RKHS of the wave functions.

VI. For the scattering of particle with spins,  the PMD-SS-optimal states are s-channel helicity conserving states (see Refs. ).

Therefore, the PMD-SS-optimal state (7) have not only the property that is the most forward-peaked quantum state but also possesses many other peculiar properties such as I-V that make it a good candidate for the description of the quantum scattering via an optimum principle. So, the systematic experimental investigation of the optimal bound (1) as well as of the bounds (34) can be of great interest since, as a direct signature of the PMD-SS-optimal state dominance in the hadron-hadron scattering phenomena, the ssaturation of the most forward-peaked limit as well as the saturation of the maximum entropy limits are expected to be experimentally observed.

4. Experimental tests

Now, in order to obtain a consistent experimental test of the optimal unitarity lower bound (1) we have compiled from the literature - the experimental data on the diffraction slope parameter $`b`$ for $`\pi ^\pm P\pi ^\pm P,`$ $`K^\pm PK^\pm P,`$ $`PPPP,`$ and $`\overline{P}P\overline{P}P`$ at all available energies.

$`\pi ^\pm P`$-scattering: In this case, for $`P_{LAB}3GeV/c`$ , the experimental data on b, $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$ are collected mainly from the original fit and data of Refs.-. To these data we added some values of $`b`$ from the linear fit of Lasinski et al.  and also calculated directly from the phase shifts analysis (PSA) of Holer et al. . In the low laboratory momenta $`P_{LAB}2`$ GeV/c, the slope parameters $`b`$ as well as the optimal slope $`b_o`$ are determined in pairs $`(b,b_o)`$, at each laboratory momenta, from the PSA . Unfortunately, the values of $`b_o`$ corresponding to the Lasinski’s data  was impossible to be calculated since the values of $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ from their original fit are not given.

$`K^\pm P`$scattering. The experimental data on b, $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$ , in the case of $`K^{}P`$, are collected from the original fit and data of Refs. , , , . For $`K^+P`$-scattering, to these data from Refs. , , ,  we added some values of $`b`$ from the linear fit of Lasinski et al.  and also those pairs $`(b,b_o)`$ calculated directly from the experimental (PSA) solutions of Arndt et al. .

$`\overline{P}P`$ and $`PP`$scatterings. The experimental data on b, $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$, in these cases, are obtained from Refs. -. The data from original fits where $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ is given numerically was first selected for to obtain the pairs $`(b,b_o)`$.

Now, for comparison of the experimental values of the slope parameter b with the optimal PMD-SS-optimal slope $`b_o`$ we used a $`(b,b_o)`$-plot of the elastic diffraction slopes as in Figs. 1a,b, where we plotted a total number $`N_p=420`$ of experimental pairs $`(b,b_o)`$. From these pairs a number $`N_p=125`$ pairs shown in Fig. 1b are coming from the experiments at $`P_{LAB}2`$ GeV/c. The dependence of the experimental slopes b and of the optimal PMD-SS-slope $`b_o`$ (see Eq. (9)) on laboratory momenta ($`P_{LAB}`$) is presented in Figs. 2-3. The values of the $`\chi ^2=_j(b_jb_{oj})^2/(ϵ_{bj}^2+ϵ_{b_oj}^2)`$, (where $`ϵ_{bj}`$ and $`ϵ_{b_oj}`$ are the experimental errors corresponding to $`b`$ and $`b_o`$ respectively) are used for the estimation of departure from the optimal PMD-SS-slope $`b_o`$, and then, we obtain the statistical parameters presented in Table 1.

5. Summary and Conclusion

The main results obtained in this paper can be summarized as follows:

(i) In this paper we proved that the optimal bound (1) is the singular solution ($`\lambda _0=0)`$ of the problem to find the unitarity lower bound on the logarithmic slope b with the constraints imposed by unitarity when $`\sigma _{el}`$ and $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ are fixed. This result is similar with that obtained recently in Ref.  for the problem to find an upper bound for the scattering entropies when $`\sigma _{el}`$ and $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ are fixed from the experimental data. So, the PMD-SS-optimal state is not only the most peaked state but also a state which saturates the maximum entropy limit;

(ii) The optimal bound (1) is verified experimentally with high accuracy (see the region b$``$b<sub>o</sub> in Fig. 1a) at all available energies for all the principal hadron-hadron scatterings.

(iii) A systematic tendency toward the saturation of the most forward-peaked optimal state limit is observed from the experimental data on the logarithmic slope of all the $`PPPP,`$ $`\overline{P}P\overline{P}P`$, $`\pi ^\pm P\pi ^\pm P`$ and $`K^\pm PK^\pm P`$ scatterings at all available laboratory momenta. The proton-proton ($`\stackrel{\mathrm{\_}}{\chi ^2}=5.06)`$ and antiproton-proton  ($`\stackrel{\mathrm{\_}}{\chi ^2}=1.86)`$ and $`K^+P`$ $`(\stackrel{\mathrm{\_}}{\chi ^2}=3.05)`$ scatterings are examples where this saturation of the most peaked optimal state limit is observed at all available laboratory momenta. In fact the validity of the principle of least distance in space of states in hadron-hadron scattering for p$`{}_{LAB}{}^{}2`$ GeV/c is well illustrated in Fig. 1b and Table 1. A similar tendency towards the saturation of the most forward peaked PMD-SS-optimal state limit is observed (see Figs. 2-3) even in the low energy regions at the energies between the resonances positions or/and in the region corresponding to the ”diffractive resonances ”.

(iv) In order to see why the experimental logarithmic slopes b$`{}_{\mathrm{exp}}{}^{}(`$ as well as the experimental entropies S<sub>L</sub> and S<sub>θ</sub>) are well described by the PMD-SS-optimal state (7) we observe that the entropy $`S_L`$ is similar to the Boltzmann entropy with a maximum value given by the logarithm of number of the optimal states. Indeed, from (31) and (34) we can conclude $`(S_LS_L^{o1}=\mathrm{ln}[(L_o+1)^2])`$ that the PMD-SS-optimal state is the state of equilibrium of the angular momenta channels considered as a quantum statistical ensemble.

One of us (D.B.I) would like to thank Prof. G. Altarelli for fruitful discussions as well as for hospitality during his stay in TH-Division at CERN

Table 1 : $`\stackrel{\mathrm{\_}}{\chi ^2}`$statistical parameters of the principal hadron-hadron scattering. In these estimations for P$`{}_{LAB}{}^{}2`$ GeV/c the errors $`ϵ_{b_o}^{PSA}(\pi ^\pm P)=0.5`$ $`GeV^2`$ and $`ϵ_{b_o}^{PSA}(K^+P)=0.3`$ $`GeV^2`$ are taken into account for the optimal slopes b<sub>o</sub> calculated from phase shifts analysis  and , respectively.

|  |  | N$`{}_{p}{}^{}(`$P$`{}_{LAB}{}^{}2GeV/c)`$ |  | N$`{}_{p}{}^{}(`$all P$`{}_{LAB}{}^{})`$ |
| --- | --- | --- | --- | --- |
| Statistical parameters | N<sub>p</sub> | $`\chi ^2/`$N<sub>p</sub> | N<sub>p</sub> | $`\chi ^2/`$N<sub>p</sub> |
| $`P`$ $`PP`$ $`P`$ | 29 | 5.01 | 32 | 5.06 |
| $`\stackrel{\mathrm{\_}}{P}P\stackrel{\mathrm{\_}}{P}P`$ | 27 | 0.56 | 45 | 1.86 |
| $`\pi ^+P\pi ^+P`$ | 12 | 2.90 | 43 | 9.70 |
| $`\pi ^{}P\pi ^{}P`$ | 14 | 2.61 | 68 | 14.9 |
| $`K^+PK^+P`$ | 20 | 0.75 | 55 | 3.05 |
| $`K^{}PK^{}P`$ | 23 | 2.18 | 97 | 6.40 |
| All experiments | 125 | 2.38 | 340 | 6.83 |

Fig. 1: The experimental values $`b`$ versus optimal values $`b_o`$ of the logarithmic slope for all principal hadron-hadron scatterings: (a) at all available momenta, and (b) only for $`P_{LAB}2`$ GeV/c. The experimental data for $`b`$,$`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$, are taken from Refs. -.

Fig. 2: The experimental values (black circles) of the logarithmic slope $`b`$ for the proton-proton and antiproton-proton scatterings are compared with the optimal PMD-SS-predictions $`b_o`$ (white circles). The experimental data for $`b`$, $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$, are taken from Refs. - (see the tex).

Fig. 3: The experimental values (black circles) of the logarithmic slope b for the principal meson-nucleon scatterings are compared with the optimal PMD-SS-predictions $`b_o`$ (white circles). The experimental data for $`b`$, $`\frac{d\sigma }{d\mathrm{\Omega }}(1)`$ and $`\sigma _{el}`$, are taken from Refs. -. The solid curves shows the values of $`b`$ calculated from the experimental PSA  (see the text).