text:
# Edwards measure and the steady state regime of a model with kinetic constraints under tapping

## 1 Introduction

The physical properties of granular materials have been extensively studied during the last few decades because of their wide presence in industry and their interest for fundamental statistical mechanics. In particular, since the compaction experiments of the Chicago group , there have been many attempts to understand the mechanism of compaction of dense powders under weak tapping and their steady state behaviour. In these systems, energy is completely dissipated after each tap and the thermal fluctuations are negligible compared to the gravitational energies involved. Such a system evolves from one blocked state to another due to the external driving without obvious detailed balance and the usual tools of statistical mechanics have to be generalized. Edwards and coworkers made the assumption that in this context it is possible to build a “thermodynamics” by using a flat measure over the blocked states in the steady state, the main ingredient being that all blocked, or metastable, configurations are equiprobable in the steady state . This is the simplest and most natural first assumption. This ergodicity in the steady state could conceivably arise from the extensive, non local nature of the tapping dynamics. It seems to be a rather hard task to define and compute the entropy of blocked states in a realistic granular medium, for example an assembly of grains in a box . Hence the Edwards measure has been recently put to the test on a wide variety of simple models. It was found in the three dimensional Kob-Andersen  and Tetris models  that the flatness assumption is good for dense systems, and also in the context of partially analytically tractable one dimensional models  and spin models on random graphs . Recent simulations on three dimensional sheared packing of spheres  have applied the Edwards measure to a realistic model and opened up the possibility of testing Edwards’ hypothesis on physical systems. Moreover, a tapping mechanism has been introduced on spin glass models and the Edwards measure was shown to be very efficient in describing phase transitions in the steady state . As the Edwards measure seemed to give very good results on the thermodynamics of the tapped Ising chain , it has also been tested on one dimensional kinetically constrained models . Variants of these models had been studied to test the Stillinger and Weber idea , which is to decompose the space of configurations into valleys, to project each valley onto its minimum (called inherent structure) and to reduce the dynamics of the system to a dynamics on inherent structures. Ising models with kinetic constraints allow one to test this decomposition in the following way : two models differ only in the constraints of the dynamics and share the same local energy minima. As the low temperature dynamics are different, the dynamics cannot be reduced to a simple sampling of inherent structures. For the same reason, the Edwards measure was expected to fail. Berg et al.  submitted these models to two kinds of tapping, which they called “thermal” and “random” and found that the Edwards measure fails, as demonstrated by De Smedt et al. in the limit of large tapping . Moreover, they argued that the validity of the flatness assumption depends on the tapping mechanism, that is to say the way energy is injected in the system. In this paper, we show how the basic Edwards measure can be generalized to build the thermodynamics of the asymptotic regime and that the validity of this measure is independent of the tapping mechanism. We shall attribute the deviation of the measure sampled during tapping simulations from the Edwards one to short distance correlations in the metastable states and try to explain why the basic flat measure is efficient at very weak tapping.

## 2 The generalized Edwards measure

Edwards’ hypothesis consists of assuming that the steady state dynamics is ergodic, i.e. the resulting measure over blocked configurations is flat. In addition, if some quantities are conserved on average, the measure must be restricted to the blocked configurations having these average quantities, as it is the case in ordinary statistical physics. Other quantities, which are not conserved, fluctuate around a value which maximizes the Edwards entropy. The original idea of Edwards and coworkers  was that an assembly of grains in a gently vibrated box is fully characterized by its density (or volume “V”), which is the only quantity to fix on average in the steady state. Hence they introduced, as a Lagrange multiplier, a “compactivity”:

$$X_{Edw}^1=\frac{S_{Edw}}{V}$$
(1)
However there is no evidence that only one quantity such as the density (or the energy for spin systems) has a non zero Lagrange multiplier and, as already mentioned , at least two quantities should be fixed on average in order to describe the steady state with a flat measure.

Let us then build a generalized Edwards measure and imagine a granular like system submitted to external forcing whose effect is to move the system from one blocked configuration to another. Let us assume that the balance between forcing and relaxation is such that in the asymptotic stationary regime exactly $`m`$ quantities $`X_k\left(k=1,\mathrm{},m\right)`$ are conserved on average. For instance, for a mixture of hard spheres of different diameters $`d_1`$ and $`d_2`$, one can consider $`X_1=h_1`$ and $`X_2=h_2`$ the mean heights of each kind of sphere. We introduce corresponding Lagrange multipliers $`\beta _k`$ and compute the grand canonical partition function:

$$Z(\{\beta _k\})=\underset{k}{}dX_ke^{_k\beta _kX_k+S(\{X_k\})}$$
(2)
where $`S(\{X_k\})`$ is the entropy of the blocked configurations restricted to that of given $`\{X_k\}`$. In the limit of a large volume or number of particles, the integrand is sharply peaked around one value $`\{X_k^{}\}`$ which maximizes $`_k\beta _kX_k+S(\{X_k\})`$. The Lagrange multipliers are given as in usual statistical mechanics by:

$$\beta _k=\frac{S}{X_k}$$
(3)
and the average of $`X_k`$ is:

$$X_k=\frac{\mathrm{log}Z}{\beta _k}=X_k^{}$$
(4)
### 2.1 The model

The model we shall consider in this paper is a variant of the Fredrickson-Anderson (FA) model, which will be refered in what follows as the symmetrically constrained Ising model (SCIM). In the original FA model , particles are deposed on a one dimensional lattice. At each site $`i`$ is associated its occupation number $`n_i=0,1`$. The total energy is $`_in_i`$ and the dynamics is constrained, that is the usual metropolis probability for a spin to flip is weighted by an acceptence ratio: $`W(n_i1n_i)=\frac{1}{2}\left(2n_{i1}n_{i+1}\right)\mathrm{min}(1,e^{\beta \mathrm{\Delta }E})`$. In this model, equilibration proceeds through elimination of isolated holes by coalescence, which is slower and slower at low temperature, as these defects are very separated, and the dynamics and the system undergoes a dynamical glass transition . This model has been studied in the context of granular compaction . It was shown to have very slow dynamics consistent with the inverse logarithmic law found in experiments , followed by a steady state well described by a flat measure over the blocked states. Here, the kinetic constraint will be changed a little: in any single move step, a particle can be added to or removed from a site only if at least one of the neighbouring sites is empty. It has been recently shown that a basic application of the Edwards measure is unable to describe the thermodynamics of the steady state of the SCIM submitted to two types of tapping :

* “random”: occupation of each site is changed with probability $`p[0,1/2]`$;
* “thermal”: one Monte-Carlo sweep is made, with Metropolis probability $`p(n_i1n_i)=(1n_{i1}n_{i+1})\mathrm{min}(1,e^{\beta \mathrm{\Delta }E})`$. This thermal tapping was introduced and studied on spin models  and lattice models  of granular matter, after the analogy between vibration and thermal noise was pointed out in .

In between taps, the system undergoes a zero temperature dynamics which corresponds to adding particles at empty sites having at least one empty neighbour, until becoming blocked in a metastable state. This dynamics can however be seen as the zero temperature Glauber dynamics of a model (without kinetic constrains) with energy per site $`E=\frac{1}{N}_i\left((1n_i)(1n_{i1}n_{i+1})n_{i1}n_{i+1}\right)`$, where only moves which strictly lower the energy are allowed. With this definition, the metastable states are now energetically metastable. The contributions of site $`i`$ to this energy is:

* $`1`$, if site $`i1`$ and $`i+1`$ are occupied;
* $`0`$, if site $`i`$ is occupied and either site $`i1`$ or site $`i+1`$ is empty;
* $`1`$ , if site $`i`$ is empty and either site $`i1`$ or site $`i+1`$ is empty.

In addition, in the following, the average occupation $`\rho =\frac{1}{N}_in_i`$, which involves no interactions, will be called “density”. This point is fundamental, as we have to keep in mind that the blocked states are reached by “gradient” descent in the energy landscape. Hence, as the basins of attraction of the metastable states are not a priori the same, the assumption that the latter are sampled in a flat manner is a strong one.

As in granular media the complexity emerges from the kinetic constraints due to hard-core repulsion and collisions, this model is thus a simple one dimensional granular medium submitted to tapping, with an energy $`E`$ driving the “falling” of the particles and a density $`\rho `$ characterizing the compacity. The definition of the entropy of metastable states is exactly that of Edwards and one can apply Edwards’ hypothesis in its original spirit.

The zero temperature dynamics stops when all empty sites are isolated. This gives a simple characterization of any metastable states as a sequence of domains of neighbouring occupied sites, separated by one empty site and allows one to predict easily the entropy or the distribution $`P(l)`$ of domain sizes among the metastable states. Our goal is then to compute as many characteristics of the steady state as possible, with a minimal set of quantities obtained by measurement, and if possible to find some circumstances where the simplest Edwards measure is a good approximation.

The description of the steady state regime by the Edwards measure fails in the regime of low density, where the average length of the domains is small. On the other hand, it seems to be fairly accurate in the high density regime, where large domains dominate. Moreover, the zero temperature dynamics involves short range interactions, so we expect that, in the blocked states sampled by the tapping dynamics, correlations at long distances are induced by correlations at short distances. In the context of the SCIM, Edwards’ hypothesis implies that all correlations are obtained from the average domain length. If one wants to improve the measure by introducing new Lagrange multipliers, one can add some which fix the average value of short length scale characteristics, such as the number of domains of length one or two.

For simplicity, and as the density $`\rho `$ and the energy per site $`E`$ are natural quantities of the model, we have computed the entropy and the distribution of domain lengths in the Edwards ensemble for given values of $`\rho `$, $`E`$ and the probability $`P(l<3)`$ that a domain has length smaller than $`3`$, which are linear combinations of $`\rho `$, $`P(l=1)`$ and $`P(l=2)`$, so short length scales are fixed on average, as explained above. To keep only two quantities, we maximize the Edwards entropy with respect to $`\alpha =P(l<3)`$ with $`\rho `$ and $`E`$ fixed. Maximizing again with respect to $`E`$ gives the simplest Edwards measure. The results are given in the appendix.

## 3 Numerical simulations

Here, we shall compare the accuracy of different generalizations of the flat measure in numerical simulations of tapping, as well as the influence of the tapping mechanism. In order to clarify how different excitations can lead to different regimes of density, it is important to separate the ingredients of the tapping (ii), that is the kinetic constraints and the thermal condition. In addition, we can combine the kinetic constraints with the random tapping (i). So we define four tapping mechanisms, depending whether, under tapping, the kinetic constraints are respected and whether the tapping is random:

* “Random Unconstrained”: occupation of each site is changed with probability $`p[0,1]`$. Notice that there is no reversal symmetry as in $`\pm 1`$ Ising spin systems, so $`p`$ can be greater than $`1/2`$;
* “Random Constrained”: one Monte-Carlo sweep is made during which the occupation of each chosen site is changed randomly with probability $`p[0,1]`$ if it has an empty neighbouring site;
* “Thermal Unconstrained”: one Monte-Carlo sweep is made with Metropolis probability $`p(n_i1n_i)=\mathrm{min}(1,e^{N\mathrm{\Delta }\rho /T_\rho })`$, $`T_\rho `$ being the tunable intensity of tapping and $`\mathrm{\Delta }\rho =\frac{12n_i}{N}`$ the variation of the density during the Metropolis step;
* “Thermal Constrained”: one Monte-Carlo sweep is made with Metropolis probability $`p(n_i1n_i)=(1n_{i1}n_{i+1})\mathrm{min}(1,e^{N\mathrm{\Delta }\rho /T_\rho })`$;

In the two latter cases, we use $`\rho `$ instead of $`E`$ in order to compare with the results of Berg et al.  (so (RU) corresponds to (i) and (TC) corresponds to (ii)) (However, this is equivalent in the (TC) case and the results are not qualitatively changed if we use $`E`$ instead of $`\rho `$ in (TU)). We have carried out simulations for each of the four tapping mechanisms above. The systems had $`N=10^5`$ and $`N=10^6`$ spins and several quantities have been recorded during $`10^6`$ taps once the steady state reached:

* the energy $`E`$ and the density $`\rho `$;
* the distribution of the domain sizes $`P(l)=Probability(\text{size}=l)`$;
* the fluctuations of $`E`$ and $`\rho `$:

$`c_E`$ $`=`$ $`N(E^2E^2)`$ (5)
$`c_\rho `$ $`=`$ $`N(\rho ^2\rho ^2)`$ (6)

As expected, $`E`$ does not maximize the entropy when $`\rho `$ only remains fixed, and a small but significant dependence of the curve of $`E`$ vs $`\rho `$ on the tapping mechanism is observed.

In order to put to the test the applicability of the canonical ensemble with three non zero Lagrange multipliers, we shall compare the distribution $`P(l)`$ of domain lengths and the fluctuations of $`\rho `$ and $`E`$ recorded during the simulations, with their corresponding values in the Edwards ensemble restricted to the configurations where energy, density and probability for domain length to be at most three, are equal to that measured during the simulations. In the following, we shall refer the corresponding measure as M3. The same procedure is carried out for two and one non zero Lagrange multipliers, where the ensemble was restricted to energy and density, or density only, with the corresponding measures referred as M2 and M1 respectively.

Let us remark that if one assumes that the distributions of the lengths of two neighbouring domains are independent, $`P(l)`$ is enough to compute all correlation functions involving a finite number of sites, so measuring $`c_\rho `$ and $`c_E`$ may be redundant. However, these fluctuations involve a large number of terms, and so are very sensitive to the deviations to the exact measure sampled during the simulation. As we shall see below, the comparisons between the different generalized measures and the numerically generated one are much more convincing when comparing the fluctuations than when comparing the distributions of domain lengths.

In figure (1) and figure (2), the computation of the fluctuations of the density obtained from the tapping simulations are displayed and compared to the ones expected from measures M1, M2 and M3, as a function of the steady-state density. As explained in the next section, different kinds of tapping cover different energy and density ranges, so we can test the generalized Edwards measure on a wide range of energies or densities. As is already known , if only $`\rho `$ is fixed, this measure is accurate only at high density. We have verified that the measure with only $`E`$ fixed works only at low energy too. Fixing both $`\rho `$ and $`E`$ gives quite good results, but there is still a difference between the tapping simulations and the value expected from the generalized Edwards measure. The distribution of domain lengths $`P(l)`$ obtained with (TC) with $`T_\rho =1.3`$ is shown in figure (3) and with (RU) with $`p=0.4`$ in figure (4). The non exponential behaviour of $`P(l)`$ at short lengths indicates that we have to fix at least two quantities on average. The computation of $`P(l)`$ using the measure M2 is better than that using M1, but a difference with the simulations remains at low density, as shown in figure (3). However, $`P(l)`$ becomes exponential as soon as $`l3`$, which indicates that the large scale degrees of liberty maximize the entropy, so that only three parameters should be enough to describe the whole distribution $`P(l)`$. So, the computation using the measure M3 is expected to predict with accuracy the fluctuations of density and the distribution of domain lengths, as it is the case in figure (1), figure (2) and figure (3). In addition, we remark that the local minimisation of energy involves three consecutive sites, so that the effective interaction length due to the kinetic constraints is three. So one can expect that if the zero temperature dynamics involves four consecutive sites, the measure M4 will be needed.

In such models, where energy is injected in the system by external forcing, there is no conservation law to insure that any given quantity must be fixed on average as a result of the equilibration between the internal relaxation into metastable states and the external driving. However, as expected, the generalized Edwards measure converges to the measure sampled during the tapping when the number of Lagrange multipliers increases. Moreover, the convergence to the original Edwards measure is more rapid as the tapping intensity is lowered and near the maximum of density, the simplest measure gives a very good approximation. At low tapping, if we consider the tapping mechanisms (TC) and (RC), which allow one to reach this regime of high density (see next section), the average domain size is large and the dynamics is dominated by the diffusion of small sequences of short domains, separated by long domains. In the language of granular media, the mobile particles are localized in regions of weak density, which are far form one another at high density and then diffuse independently. This is reminiscent of similar results in the context of the Kob-Andersen model, for which Edwards’ hypothesis was found to apply at high density. In these models, as well as in granular media, in the very compact regime, the majority of particles are unable to move during the taps because of the hard core constraints.

Hence, we conjecture that this scenario is more general: let us consider a granular like system, that is an assembly of hard “heavy” particles, evolving among blocked configurations thanks to a macroscopic forcing. If the driving excitation is weak enough so that in the steady state the defects (regions where the density is low) are distant one from each others, the dynamical measure over the blocked configurations is flat. On the contrary, if the system is near the random loose packing, it is very heterogeneous in space and a majority of particles are allowed to move, contributing to large avalanches which break the ergodicity.

## 4 Comparing different tapping mechanisms

It has been argued through the tapping mechanisms (i) and (ii) that “thermal” tapping is much more efficient in sampling the configurations in the flat manner than the “random” tapping. Moreover, it was added that the former allows the system to reach high densities, whereas the latter was confined below $`\rho ^{}=3/4`$. Here we shall clarify this issue by separating the influence of the “thermal” or “random” nature of the excitation and the presence or absence of a kinetic constraint during the tap by comparing the results of the simulations with the mechanisms (TC), (RC), (TU) and (RU).

Indeed, if the same kinetic constraint as that of the zero temperature dynamics is not imposed during the tap, domains can split or coalesce as the system is excited, whereas the number of domains changes during the relaxation only through nucleation. Thus large domains are unstable with respect to (TU) and (RU) and stable with respect to (TC) and (RC). This remark allows one to compute the maximal steady state density accessible to (TU) and (RU), obtained in the limit of zero excitation intensity. To do so, we assume that the tapping is so weak that in a given sequence of sites, only one site is changed. Neglecting the correlations of the lengths of consecutive domains, we focus on three consecutive domains, where at most one change occurs in the central one (which size is $`l=_llP(l)`$) or at its frontier during a tap. The average density is given by $`\rho =l/(l+1)`$ and its variation after one tap (e.g. in the next metastable state) is:

$$N\mathrm{\Delta }\rho =\left((N\rho +1)p(a)+N\rho p(b)+(N\rho 1)p(c)\right)N\rho $$
(7)
where

$`p(a)`$ $`=`$ $`{\displaystyle \frac{1}{l+1}}`$ (8)
$`p(b)`$ $`=`$ $`{\displaystyle \frac{2}{l+1}}`$
$`p(c)`$ $`=`$ $`{\displaystyle \frac{l2}{l+1}}`$
and $`p(a)`$,$`p(b)`$ and $`p(c)`$ are the probability of the contributions to the variation of the density displayed on figure (5). Gathering these three terms gives:

$$N\mathrm{\Delta }\rho =\frac{3l}{1+l}$$
(9)
Here, because any empty site is shared by two domains in a metastable state, we do not add a particle in the empty site at the left of the central domain to avoid redundancies and then the denominator in Eq. (9) is just the number of possible moves. This gives $`\rho ^{}=3/4`$ in the large time limit, even if the system is prepared in a high density state.

On the contrary, if the dynamics is constrained, the zero tapping limit corresponds to moving domain walls and then slowly eliminating small domains as in the zero temperature evolution of the one dimensional Ising model . Hence, the mean length of the domain walls growths until it is of the order of the size of the system, and the density approaches very slowly the maximum value possible.

Now, it is clear from figure (1) and figure (2) that (TU) and (RU) on one hand, (TC) and (RC) on the other hand lead to comparable deviations of the the measures M1 and M2 from the dynamical one. So, the validity of Edwards’ hypothesis is independent of the way the system is tapped. This is not surprising, since all configurations are connected by the unconstrained part of the four tapping mechanisms (thermal or random) so that the ergodicity is broken by the kinetic constraint and the zero temperature dynamics. Hence, the only relevant difference between the four tapping mechanism introduced here is whether or not the kinetic constraint is respected during the excitation. However, the only influence of this constraint on the accuracy of the Edwards measure is through the range of densities accessible.

Even if the measures M1, M2, or M3 are insensitive to the kind of excitation, we can use the results of the tapping simulations with mechanisms (TC) and (RC), or (TU) and (RU) to find whether the measure M1 is accurate or not, without knowing a priori the Edwards entropy. Indeed, as far as lattice models as the SCIM involved here, the Edwards measure M1 can be computed at least numerically and compared to one obtained dynamically. However, as far as realistic granular media are concerned, this is no longer possible. If we compare the values of the observables measured during the tapping simulations with (TC) and (RC), the fluctuations of density for instance, we find some small differences, in the regime of density where M1 does not apply. One can explain these differences by considering the measure M2. Indeed, if for instance the density and the energy are fixed, the value of the entropy in the steady state, when the tapping intensity decreases, is a path on a two dimensional surface, which depends on the tapping mechanism: as the energy is injected in the system in two different manners, the Lagrange multipliers are not expected to be the same for the same value of the density. On the contrary, as M1 involves only the density, characteristics of the steady state, such as density fluctuations, should not depend on the excitation mechanism for a given value of the steady state density, if M1 applies. So one can try to imagine how the Edwards measure can be put to the test experimentally. Let us assume that a given steady state packing fraction $`\varphi `$ of the same grains can be obtained by several forcing mechanisms, like shear and vibration for instance. If the amplitude of the fluctuations of the packing fraction differ for different kinds of excitation, $`\varphi `$ cannot be the only relevant macroscopic quantity.

## 5 Conclusion

In this paper, we have addressed the possibility of describing the steady state regime of a simple model of granular matter by using a flat measure. In the high density regime, the knowledge of the mean length of domains was enough to give a qualitative description of global quantities, like the fluctuations of $`\rho `$. However, Lagrange multipliers corresponding to shortest length scales have been introduced in order to compute with accuracy all the quantities of interest for all the densities accessible. We have shown that by introducing a small number (related to the length scale of the effective interaction involved in the zero temperature dynamics) of multipliers, one may characterize the system in cases where the basic flat measure fails. Moreover, we have given a more general context where the latter applies. Its success in the weak tapping limit was attributed to the diffusion of the regions where the density is low, well separated in space. We then compared different tapping mechanisms and showed that the Edwards measure was indifferent to the way the system was excited. In addition, we have proposed a way of testing the original Edwards measure without any a priori information about the entropy of metastable states, by applying different kinds of tapping.

I am grateful to D. S. Dean and A. Barrat for numerous corrections in the manuscript and stimulating discussions.

## Appendix

The entropy $`s`$ and the distribution $`P(l)`$ of domain sizes in different generalized Edwards ensembles can be obtained by simple combinatorial arguments. Especially, if one imposes that a given domain has size $`l`$ and counts the number of possible metastable states which fulfill this constraint (with $`\rho `$ or $`E`$ or $`\alpha `$ fixed), the leading term gives $`e^{Ns}`$ and the term of order $`1`$ gives $`P(l)`$.

### A.1 The original measure (M1)

The entropy of metastable states per site is:

$$s(\rho )=\rho \mathrm{log}\frac{2\rho 1}{\rho }+(1\rho )\mathrm{log}\frac{2\rho 1}{1\rho }$$
(10)
and the distribution of domain lengths is exponential:

$$P(l)=\frac{1\rho }{\rho }\left(\frac{2\rho 1}{\rho }\right)^{l1}\theta (l1)$$
(11)
where $`\theta `$ is the Heaviside step function.

### A.2 Fixing two quantities (M2)

The entropy of metastable states per site is:

$`s(\rho ,E)`$ $`=`$ $`(1\rho )\mathrm{log}(1\rho )+(2\rho 1)\mathrm{log}(2\rho 1)2(\rho +E)\mathrm{log}(\rho +E)`$
$``$ $`(12\rho E)\mathrm{log}(12\rho E)(\rho 1E)\mathrm{log}(\rho 1E)`$
and the distribution of domain lengths is exponential only for $`l2`$:

$$P(l)=\left(\frac{12\rho E}{\rho +E}\right)^{\theta (l2)}\frac{12\rho E}{1\rho }\frac{\rho 1E}{\rho +E}\left(\frac{\rho 1E}{2\rho 1}\right)^l\theta (l1)$$
(13)
If one maximizes the entropy with respect to $`E`$, then $`E=\frac{(12\rho )^2}{\rho }`$ and the calculations with only $`\rho `$ fixed is recovered.

### A.3 Fixing three quantities(M3)

The entropy of metastable states per site is

$`s(\rho ,E,\alpha )`$ $`=`$ $`(\rho 1E)\mathrm{log}(\rho 1E)+(1\rho )\mathrm{log}(1\rho )`$
$``$ $`(\rho 2E+\alpha )\mathrm{log}(\rho 2E+\alpha )2(1\alpha )\mathrm{log}(1\alpha )`$
$``$ $`(\rho +E1+\alpha )\mathrm{log}(\rho +E1+\alpha )`$
and the distribution of domain lengths is exponential only for $`l3`$:

$`P(l)`$ $`=`$ $`{\displaystyle \frac{\rho +E1+\alpha }{(1\alpha )(12\rho E)}}\left[{\displaystyle \frac{(\rho 1E)(\rho +E1+\alpha )}{(\rho 2E+\alpha )(12\rho E)}}\right]^{\theta (l2)}`$
$`\times `$ $`\left[{\displaystyle \frac{(1\alpha )^2}{(\rho +E1+\alpha )(\rho 2E+\alpha )}}\right]^{\theta (l3)}\left[{\displaystyle \frac{\rho 2E+\alpha }{\rho 1E}}\right]^l\theta (l1)`$
Here again, the value of $`\alpha `$ which maximizes the entropy is $`\alpha =1\frac{(\rho 1E)(\rho +E)}{2\rho 1}`$ and with this value the calculations with $`\rho `$ and $`E`$ fixed are recovered.

## References


format: text

text:
# Hydrodynamic Model for Particle Size Segregation in Granular Media

## I Introduction

The physics of granular materials is a subject of current interest. A granular medium is a system of many macroscopic heterogeneous particles with dissipative interactions. One of the outstanding problems is the so–called “Brazil Nut effect”: When a large intruder particle placed at the bottom of a vibrated bed tends to the top. This size segregation is due to the nonequilibrium, dissipative nature of granular media. Granular materials are handled in many industries. Many industrial machines that transport granular materials use vertical vibration to fluidize the material, and the quality of many products is affected by segregation. Size segregation is one of the most intriguing phenomena found in granular physics. A deeper understanding of this effect is therefore interesting for practical applications, and also represent theoretical challenge.

A series of experiments and computer simulations have elucidated different size segregation mechanisms, including vibration frequency and amplitude; particle size and size distribution; particle shape; and other properties such as density and elastic modulus.

Several possible mechanisms for size segregation have been proposed. One is segregation in presence of convection observed experimentally in three dimensions by Knight et al.,, and by Duran al., in two dimensions, under conditions of low amplitude and high acceleration vibration. In this case, both intruder and the small particles are driven up along the middle of the cell, and while the smaller particles are carried down in a convection roll near the walls the intruder remains trapped on the top. In experiments performed by Vanel et al., in three dimensions, they observed two convective regimes separated by a critical frequency. The first regime is associated with heaping and the second regime is similar to the one observed in Ref.. Also, they reported a nonconvective regime observing a size dependent rise velocity. Employing large molecular dynamics simulations in two dimensions Pöschel and Herrmann, and in three dimension Gallas et al.,, have recovered several aspects that are seen in experiments and recognize the lack of a theoretical description of the exact mechanism driving the segregation and the role of convection.

Other segregation mechanism is associated to the percolation of small grains. Based on a Monte Carlo computer simulation, Rosato et al., argue that each cycle of the applied vibration causes all the grains to detach from the base of the container. Then, the smaller particles fall relatively freely, while the larger particles require larger voids to fall downwards. The large grains therefore effectively rise through the bed. In the context of large–amplitude, low–frequency vertical shaking process (tapping), Jullien et al., predict a critical size ratio below which segregation does not occur. This provoked some controversy and this threshold may be an artifact of the simulation model based on the “steepest descent algorithm”. Experiments in Helle Shaw cells observed an intruder size dependent behavior, where the segregation rate increases with the size ratio between the intruder and the surrounding particles. Duran et al., formulated a geometrical theory for segregation based on the arching effect. They also claim experimental evidence for a segregation size threshold. In this picture the intruder contributes to the formation of an arch sustained on small grains on both sides. Between each agitation the small particles tend to fill the region below the arch. So, at each cycle the small particles move downward and the intruder effectively rises. Using a modification of the algorithm proposed by Rosato et al., Dippel and Luding find a good qualitative agreement with the non–convective and size–dependent rising.

In another context, Caglioti et al., considered the geometrical properties of mixtures in the presence of compaction. They established a relation for the effective mobilities of different particles in heterogeneous situations.

The effect of the intruder density was studied by Shinbrot and Muzzio. They observed an oscillating motion of the intruder on the top, which corresponds to the “whale effect” predicted by Pöschel and Herrmann. Also, they observed a reverse buoyancy in shaken granular beds. Ohtsuki et al., performed molecular dynamics simulations in two dimensions and studied the effects of intruder size and density on the height, and found no segregation threshold. Recently Shishoda and Wassgren performed two dimensional simulations to model segregation in vibrofluidized beds. They reported an height dependence with the density ratio between the intruder and the surrounding particles. In their model the intruder position results from a balance between the granular pressure (buoyant force) within the bed and the intruder weight. Their approach is in some sense similar to the model that we propose in this article.

Subject to an external force, granular materials locally perform random motions as a result of collisions between grains, much like the molecules in a gas. This picture has inspired several authors to use kinetic theories to derive continuum equations for the granular flow–field variables. Some of these theories have been generalized to multicomponent mixtures of grains. For different size particles in the presence of a temperature gradient, Arnarson and Willits, found that larger, denser particles tend to be more concentrated in cooler regions. This result was confirmed by numerical simulations. However, this mechanism of segregation is a natural consequence of the imposed gradient of temperature and its not related to the nature of the grains.

In this article we address the problem of size segregation using a kinetic theory approach in two and three dimensions ($`D=2,3`$). We consider the case of an intruder particle immersed in a granular bed. We propose a segregation mechanism based on the difference of densities between different regions of the system, which give origin to a buoyant force that acts on the intruder. The difference of densities is caused by the difference between the mean kinetic energy among the region around the intruder and the medium without intruder. The dissipative nature of the collisions between the particles of a granular media is responsible for this mean energy difference, and modifies the mobility of the particles.

The plan of this article is as follows. In Sec. II we derive a continuum formulation for the granular fluid, and introduce the definition of the “granular temperature”. In Sec. III we propose an analytic method to estimate the local temperature in the system. In Sec. IV we introduce the coefficient of thermal expansion. In Sec. V explicit solutions of the time dependence of height and velocity of the large particle are calculated. We can explain the size ratio dependence of the rise velocity and address the issue of the critical size ratio to segregation. To validate our arguments we make comparisons with previous experimental data.

## II Continuum formulation

We consider an intruder particle of mass $`m_I`$ and radius $`r_I`$ immersed in a granular bed. The granular bed is formed of $`N`$ monodisperse particles of mass $`m_F`$ and radius $`r_F`$. The particles are modeled by inelastic hard disks ($`D=2`$) or spheres ($`D=3`$) in a $`D`$–dimensional volume $`V=L^D`$ of size $`L`$. The size ratio is denoted $`\varphi =r_I/r_F`$. The particles interact via binary encounters. The inelasticity is specified by a restitution coefficient $`e1`$. We assume this restitution coefficient to be a constant, independent on the impact velocity and the same for the fluid particles and the intruder. The post collisional velocities $`𝐯^{}`$ are given in terms of the pre–collisional velocities $`𝐯`$ by

$$𝐯_{1,2}^{}=𝐯_{1,2}\frac{m_{red}(1+e)}{m_{1,2}}[(𝐯_1𝐯_2)\widehat{𝐧}]\widehat{𝐧},$$
(1)
where the labels $`1`$ or $`2`$ specify the particle, $`\widehat{𝐧}`$ is the unit vector normal to the tangential contact plane pointing from $`1`$ to $`2`$ at the contact time, and the reduced mass $`m_{red}=m_1m_2/(m_1+m_2)`$. To calculate the dissipated energy we consider that energy is dissipated only by collisions between pairs of grains. In a binary collision the energy dissipated is proportional to $`\mathrm{\Delta }E=m_{red}(1e^2)v^2/2`$, where $`v`$ is the mean velocity of the particles.

In this work we use a generalized notion of temperature. In a vibrofluidized granular material a “granular temperature” $`T_g`$ can be defined to describe the random motion of the grains and is the responsible for the pressure, and the transport of momentum and energy in the system. The granular temperature $`T_g`$ is defined proportional to the mean kinetic energy $`E`$ associated to the velocity of each particle

$$\frac{D}{2}T_g=\frac{E}{N}=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(\frac{1}{2}m_iv_i^2\right).$$
(2)
We expect a continuum limit to hold for $`N1`$, when the small particles may be considered as forming a granular fluid. In order to develop an analytic study, we assume that the uniformly heated granular fluid can be described by the standard hydrodynamic equations derived from kinetic theories for granular systems. In this study, we focus on a steady state with no macroscopic flow.

The balance equation for the energy is

$$𝐪=\gamma ,$$
(3)
where $`𝐪`$ is the flux of energy and $`\gamma `$ is the average rate of dissipated energy due to the inelastic nature of the particles collisions. The constitutive relation for the flux of energy,

$$𝐪=\kappa T_g,$$
(4)
defines the thermal conductivity $`\kappa `$. Consequently, we have

$$\left(\kappa T_g\right)=\gamma .$$
(5)
A uniformly fluidized state can be realized when the granular material is vibrated in the vertical direction, typically as $`z_0(t)=A_0\mathrm{sin}(\omega _0t)`$, with the amplitude $`A_0`$ and the frequency $`\omega _0=2\pi f`$, so that one can define a typical velocity $`u_0=A_0\omega _0`$. In the experiments the excitation is described by the dimensionless amplitude $`\mathrm{\Gamma }_0=A_0\omega _0/g`$, where $`g`$ is the gravitational acceleration. As a first approximation the effect of the external force experienced by the particles due to the gravitational field is neglected in the description of the granular flow. Experimentally this corresponds to the regime $`\mathrm{\Gamma }_01`$. So, the momentum balance, in the steady state, implies that the pressure $`p`$ is constant throughout the system.

The hydrodynamic equations close with the state equation, the collisional dissipation $`\gamma `$ and the transport coefficients for a granular medium. In the limit $`N1`$ the constitutive relations are determined as function of the properties of the small grains. The transport coefficients are assumed to be given by the Enskog theory for dense gases in the limit of small inelasticity.

The total pressure should be essentially equal to that of the small particles, the contribution of the intruder being negligible, since $`N1`$. For a dense system the pressure is related to the density by the virial equation of state, which in the case of inelastic particles is

$$p=\frac{1+e}{2}nT_g\left[1+\frac{\mathrm{\Omega }_D}{2D}ng_0(2r_F)^D\right],$$
(6)
where $`n=N/V`$ is the number density of small grains, $`\mathrm{\Omega }_D=2\pi ^{D/2}/\mathrm{\Gamma }(D/2)`$ is the surface area of a $`D`$–dimensional unit sphere, $`\nu `$ is the area (volume) fraction $`\nu =\mathrm{\Omega }_Dnr_F^D/D`$, and $`g_0`$ is the pair correlation function for two fluid particles. In two dimensions the pair correlation function is

$$g_0=\frac{\left(1\frac{7}{16}\nu \right)}{\left(1\nu \right)^2},$$
(7)
with the area fraction $`\nu =n\pi r_F^2`$. In three dimensions the pair correlation function is

$$g_0=\frac{\left(2\nu \right)}{2\left(1\nu \right)^3},$$
(8)
with the volume fraction $`\nu =4\pi nr_F^3/3`$.

The state–dependent thermal conductivity possesses the general form

$$\kappa =\kappa _0\sqrt{T_g},$$
(9)
where the prefactor $`\kappa _0`$ is a function of the fluid particle properties, and can be calculated using a Chapman–Enskog procedure through the solution of Enskog transport equation. The explicit expressions of these prefactors are given in Appendix A.

To estimate the collisional dissipation rate $`\gamma `$ we consider the loss of average kinetic energy per collision and per unit time. In a binary collision the kinetic energy dissipated can be expressed in terms of the granular temperature as $`\mathrm{\Delta }T_g=(1e^2)T_g/2`$. For the fluid particles, the average collision frequency $`\omega _F`$ is proportional to $`\omega _F\sqrt{T_g}`$, and we assume that it is given by the Enskog collision frequency

$$\omega _F=\frac{\mathrm{\Omega }_D}{\sqrt{2\pi }}ng_0(2r_F)^{D1}\left(\frac{2}{m_F}\right)^{1/2}T_g^{1/2}.$$
(10)
This form for the frequency of collisions is justified for a granular medium. This is a consequence that the average spacing between nearest neighbor $`s`$ is supposed to be less than the grain diameter $`(s2r_F)`$. Multiplying $`\mathrm{\Delta }T_g`$ by the collision rate $`\omega _F`$ and the number density $`n=N/V`$, we obtain the collisional dissipation rate $`\gamma _F`$ for the fluid particles

$$\gamma _F=\frac{\mathrm{\Omega }_D}{2\sqrt{2\pi }}(1e^2)n^2g_0(2r_F)^{D1}\left(\frac{2}{m_F}\right)^{1/2}T_g^{3/2}.$$
(11)
In order to simplify the mathematical notation let us express $`\gamma _F`$ as

$$\gamma _F=\xi _FT_g^{3/2},$$
(12)
where the dissipation factor $`\xi _F`$ contains the prefactors which multiply $`T_g^{3/2}`$ in Eq.(11), this is

$$\xi _F\frac{\mathrm{\Omega }_D}{2\sqrt{2\pi }}(1e^2)n^2g_0(2r_F)^{D1}\left(\frac{2}{m_F}\right)^{1/2}.$$
(13)
To understand the essential features of the intruder’s presence in the granular medium, it is adequate to adopt a simplified point of view. If the mean velocity of the fluid particles is $`u`$, the flux of fluid particles which strikes the intruder’s surface can be estimate as $`nu`$. Multiplying this flux by the area of the intruder $`\mathrm{\Omega }_Dr_I^{(D1)}`$, we can calculate the number of fluid particles which strike the surface of the intruder per unit time, and written in terms of the granular temperature we have

$$\omega _I=\frac{\mathrm{\Omega }_D}{\sqrt{2\pi }}ng_0r_I^{D1}\left(\frac{m_I+m_F}{m_Im_F}\right)^{1/2}T_g^{1/2}.$$
(14)
So, the local density of kinetic energy dissipated in the region near the intruder is

$$\gamma _I=\frac{\mathrm{\Omega }_D}{2\sqrt{2\pi }}(1e^2)\frac{n}{V}g_0(r_F+r_I)^{D1}\left(\frac{m_I+m_F}{m_Im_F}\right)^{1/2}T_g^{3/2}.$$
(15)
In the simplified form Eq.(15) can be expressed as

$$\gamma _I=\xi _IT_g^{3/2},$$
(16)
where the dissipation factor $`\xi _I`$ is defined as

$$\xi _I\frac{\mathrm{\Omega }_D}{2\sqrt{2\pi }}(1e^2)\frac{n}{V}g_0(r_F+r_I)^{D1}\left(\frac{m_I+m_F}{m_Im_F}\right)^{1/2}.$$
(17)
## III Local temperature difference

The intruder’s presence modifies the local temperature of the system due to the collisions that happen at its surface. The number of collisions on the surface increases with the size of the particle, but the local density of dissipated energy diminishes. From Eq.(5) we can calculate within a sphere of radius $`r_0`$ the value of the temperature in the granular fluid in presence of the intruder and compare it with the temperature in the granular fluid without intruder, we will denote these temperatures $`T_1`$ and $`T_2`$ respectively (see Fig.(1)). This is a simple method to estimate the temperature difference between a region with intruder and a region without intruder $`\mathrm{\Delta }T_g=T_1T_2`$.

Let us concentrate on solutions with radial symmetry. The solutions of Eq.(5), for an arbitrary dimension $`D`$, satisfy the equation

$$\frac{1}{r^{D1}}\frac{d}{dr}\left(r^{D1}\kappa _0T_g^{1/2}\frac{dT_g}{dr}\right)=\xi T_g^{3/2}.$$
(18)
This nonlinear differential equation can be simplified by the fact that the pressure is considered constant throughout the system and remembering that $`pT_g`$. So, linearizing Eq.(18) the resulting equation may be written in terms of $`wT_g^{1/2}`$,

$$\frac{d^2w}{dr^2}+\frac{D1}{r}\frac{dw}{dr}=\lambda ^2w,$$
(19)
where

$$\lambda ^2\frac{\xi }{2\kappa _0}.$$
(20)
The collisional dissipation rate can be decomposed in two parts. We propose this decomposition supposing that the energy dissipation around the intruder is dominated by the collisions between the small grains and the intruder, then the dissipation rate in this region is given by Eq.(15). In the rest of the system the dissipation rate is dominated by the collisions between small grains only. In this case the dissipation is given by Eq.(11).

First, let us consider the “inhomogeneous case” when the intruder is localized in the center of the system ($`r=0`$), see Fig.(1a). The dissipation factor $`\xi `$ can be decomposed in two parts: $`\xi =\xi _I`$ for the region near the intruder $`(r=r_I)`$, and $`\xi =\xi _F`$ for the region $`(r_I<rr_0)`$, where $`r_0`$ is the radius of the considered region.

For the inhomogeneous case we express Eq.(19) as

$$\frac{d^2w}{dr^2}+\frac{D1}{r}\frac{dw}{dr}=\{\begin{array}{cc}\lambda _I^2w\hfill & \text{for }0<rr_I,\hfill \\ & \\ \lambda _F^2w\hfill & \text{for }r_I<rr_0,\hfill \end{array}$$
(21)
where, $`\lambda _I^2\xi _I/2\kappa _0`$ and $`\lambda _F^2\xi _F/2\kappa _0`$. The solution of Eq.(21) is determined by the boundary conditions imposed upon the system. As boundary condition we suppose that the system is enclosed by an external surface of radius $`r_0`$ at temperature $`T_g(r_0)=T_0`$ (respectively, $`w(r_0)=w_0`$).

Let us denote $`T_1(r)`$ the granular temperature for the region $`(0<rr_I)`$, and $`T_{1+}(r)`$ the granular temperature for the region $`(r_I<rr_0)`$ (respectively, $`w_{}(r)\sqrt{T_1(r)}`$ and $`w_+(r)\sqrt{T_{1+}(r)}`$). The intruder’s presence imposes internal boundary conditions. On the inner surface, the temperature should satisfy

$$w_{}(r)|_{r=r_I}=w_+(r)|_{r=r_I}.$$
(22)
The flux of energy also imposes another internal boundary condition. If we suppose the flux of energy continuous on the inner surface, from Eq.(4) the granular temperature should satisfy

$$\frac{dw_{}(r)}{dr}|_{r=r_I}=\frac{dw_+(r)}{dr}|_{r=r_I}.$$
(23)
### A Solution for 2D

The solutions to Eq.(21) for $`D=2`$ are a linear combination of the modified Bessel function of order zero $`w_1(r)=\{I_0(\lambda r),K_0(\lambda r)\}`$. The general solution is

$$w_{}(r)=A_{}I_0(\lambda _Ir)+B_{}K_0(\lambda _Ir)\text{ for }0<rr_I,$$
(24)
and

$$w_+(r)=A_+I_0(\lambda _Fr)+B_+K_0(\lambda _Fr)\text{ for }r_I<rr_0,$$
(25)
where $`A_{}`$, $`A_+`$, $`B_{}`$ and $`B_+`$ are constants that must be determined from the boundary conditions.

The function $`K_0(\lambda r)`$ diverges when $`r0`$, then

$$B_{}=0.$$
(26)
When $`r=r_0`$ the Eq.(25) should satisfy the boundary condition

$$w_+(r)|_{r=r_0}=w_0,$$
(27)
this is,

$$A_+I_0(\lambda _Fr_0)+B_+K_0(\lambda _Fr_0)=w_0.$$
(28)
On the inner surface the boundary condition (22) $`w_{}(r_I)=w_+(r_I)`$ leads to

$`A_{}I_0(\lambda _Ir_I)`$ $`=`$ $`A_+I_0(\lambda _Fr_I)+B_+K_0(\lambda _Fr_I)`$ (29)
$`A_{}`$ $`=`$ $`A_+{\displaystyle \frac{I_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}}+B_+{\displaystyle \frac{K_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}}.`$ (30)
The inner boundary condition (23) leads to

$`A_{}\lambda _II_1(\lambda _Ir_I)`$ $`=`$ $`A_+\lambda _FI_1(\lambda _Fr_I)B_+\lambda _FK_1(\lambda _Fr_I)`$ (31)
$`A_{}`$ $`=`$ $`\left({\displaystyle \frac{\lambda _F}{\lambda _I}}\right)\left[A_+{\displaystyle \frac{I_1(\lambda _Fr_I)}{I_1(\lambda _Ir_I)}}B_+{\displaystyle \frac{K_1(\lambda _Fr_I)}{I_1(\lambda _Fr_I)}}\right].`$ (32)
Equating Eqs.(30) and (32) we find

$`{\displaystyle \frac{A_+}{B_+}}`$ $`=`$ $`{\displaystyle \frac{\lambda _FI_0(\lambda _Ir_I)K_1(\lambda _Fr_I)+\lambda _II_1(\lambda _Ir_I)K_0(\lambda _Fr_I)}{\lambda _FI_0(\lambda _Ir_I)I_1(\lambda _Fr_I)\lambda _II_1(\lambda _Ir_I)I_0(\lambda _Fr_I)}},`$ (33)
$``$ $`\mathrm{\Theta }_{AB}`$ (34)
From Eqs. (28) and (34) the constant $`B_+`$ should be

$$B_+=\frac{w_0}{\mathrm{\Theta }_{AB}I_0(\lambda _Fr_0)+K_0(\lambda _Fr_0)}.$$
(35)
Substituting Eq.(35) into (34) we have

$$A_+=\frac{w_0\mathrm{\Theta }_{AB}}{\mathrm{\Theta }_{AB}I_0(\lambda _Fr_0)+K_0(\lambda _Fr_0)}.$$
(36)
Substituting Eqs.(35) and (36) into (30) we have

$$A_{}=\frac{w_0}{\mathrm{\Theta }_{AB}I_0(\lambda _Fr_0)+K_0(\lambda _Fr_0)}\left[\mathrm{\Theta }_{AB}\frac{I_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}+\frac{K_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}\right].$$
(37)
The granular temperature in the inhomogeneous case is

$$T_1(r)=\{\begin{array}{cc}(A_{}I_0(\lambda _Ir))^2\hfill & \text{for }0<rr_I,\hfill \\ & \\ (A_+I_0(\lambda _Fr)+B_+K_0(\lambda _Fr))^2\hfill & \text{for }r_I<rr_0,\hfill \end{array}$$
(38)
where the constant $`A_{}`$, $`A_+`$ and $`B_+`$ are given by the Eqs.(37), (36) and (35), respectively.

In the “homogeneous case”, see Fig.(1b), the prefactor $`\lambda _I=0`$. Then the granular temperature $`T_2(r)`$ is

$$T_2(r)=\left(\frac{I_0(\lambda _Fr)}{I_0(\lambda _Fr_0)}\right)^2T_0.$$
(39)
Now we are interested in determining the temperature difference $`\mathrm{\Delta }T_g`$ between case $`1`$ and $`2`$ in the granular fluid. For this we calculate the granular temperatures at $`r=0`$. When $`r0`$ the modified Bessel function of zero order tends to $`1`$. So, in Eq.(38) and (39) in the limit $`r0`$ we find that

$$\begin{array}{cc}(I_0(\lambda _Ir))^21,\hfill & \\ (I_0(\lambda _Fr))^21.\hfill & \end{array}$$
(40)
Then, the temperature difference is

$$\mathrm{\Delta }T_g=\left[\left(\frac{1}{\mathrm{\Theta }_{AB}I_0(\lambda _Fr_0)+K_0(\lambda _Fr_0)}\left[\mathrm{\Theta }_{AB}\frac{I_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}+\frac{K_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}\right]\right)^2\left(\frac{1}{I_0(\lambda _Fr_0)}\right)^2\right]T_0$$
(41)
in two dimensions.

### B Solution for 3D

When $`D=3`$, the solution of Eq.(21) is given in terms of the spherical modified Bessel functions of zero order $`w_1(r)=\{i_0(\lambda r)=\mathrm{sinh}(\lambda r)/\lambda r,k_0(\lambda r)=e^{\lambda r}/\lambda r\}`$. The general solution in this case is

$$w_{}(r)=A_{}i_0(\lambda _Ir)+B_{}k_0(\lambda _Ir)\text{ for }0<rr_I,$$
(42)
and

$$w_+(r)=A_+i_0(\lambda _Fr)+B_+k_0(\lambda _Fr)\text{ for }r_I<rr_0.$$
(43)
The function $`k_0(\lambda r)`$ diverges when $`r0`$, then

$$B_{}=0.$$
(44)
The constants $`A_{}`$, $`A_+`$ and $`B_+`$, are calculated from the boundary conditions in a similar way as before.

$`A_{}`$ $`=`$ $`{\displaystyle \frac{w_0}{\mathrm{\Theta }_{AB}i_0(\lambda _Fr_0)+k_0(\lambda _Fr_0)}}\left[\mathrm{\Theta }_{AB}{\displaystyle \frac{i_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}}+{\displaystyle \frac{k_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}}\right],`$ (45)
$`A_+`$ $`=`$ $`{\displaystyle \frac{w_0\mathrm{\Theta }_{AB}}{\mathrm{\Theta }_{AB}i_0(\lambda _Fr_0)+k_0(\lambda _Fr_0)}},`$ (46)
$`B_+`$ $`=`$ $`{\displaystyle \frac{w_0}{\mathrm{\Theta }_{AB}i_0(\lambda _Fr_0)+k_0(\lambda _Fr_0)}},`$ (47)
where in this case the factor $`\mathrm{\Theta }_{AB}`$ is

$$\mathrm{\Theta }_{AB}=\frac{\lambda _Fi_0(\lambda _Ir_I)k_1(\lambda _Fr_I)+\lambda _Ii_1(\lambda _Ir_I)k_0(\lambda _Fr_I)}{\lambda _Fi_0(\lambda _Ir_I)i_1(\lambda _Fr_I)\lambda _Ii_1(\lambda _Ir_I)i_0(\lambda _Fr_I)}.$$
(48)
The granular temperature in the inhomogeneous case in 3D is

$$T_1(r)=\{\begin{array}{cc}(A_{}i_0(\lambda _Ir))^2\hfill & \text{for }0<rr_I,\hfill \\ & \\ (A_+i_0(\lambda _Fr)+B_+k_0(\lambda _Fr))^2\hfill & \text{for }r_I<rr_0,\hfill \end{array}$$
(49)
where the constant $`A_{}`$, $`A_+`$ and $`B_+`$ are given by the Eqs.(45), (46) and (47).

In the “homogeneous case” the prefactor $`\lambda _I=0`$. Then the granular temperature $`T_2(r)`$ is

$$T_2(r)=\left(\frac{i_0(\lambda _Fr)}{i_0(\lambda _Fr_0)}\right)^2T_0.$$
(50)
Again the temperature difference $`\mathrm{\Delta }T_g`$ is calculate at $`r=0`$ between case $`1`$ and $`2`$. When $`r0`$ the spherical Bessel function of zero order tends to $`1`$, then

$$\begin{array}{cc}\left(\frac{\mathrm{sinh}(\lambda _Ir)}{\lambda _Ir}\right)^21,\hfill & \\ \left(\frac{\mathrm{sinh}(\lambda _Fr)}{\lambda _Fr}\right)^21.\hfill & \end{array}$$
(51)
Then, the temperature difference is

$$\mathrm{\Delta }T_g=\left[\left(\frac{1}{\mathrm{\Theta }_{AB}i_0(\lambda _Fr_0)+k_0(\lambda _Fr_0)}\left[\mathrm{\Theta }_{AB}\frac{i_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}+\frac{k_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}\right]\right)^2\left(\frac{1}{i_0(\lambda _Fr_0)}\right)^2\right]T_0$$
(52)
in three dimensions.

### C Energy equipartition breakdown

Let us define the temperature ratio $`\tau T_1(0)/T_2(0)`$. In two dimension we have

$$\tau =\left(\frac{I_0(\lambda _Fr_0)}{\mathrm{\Theta }_{AB}I_0(\lambda _Fr_0)+K_0(\lambda _Fr_0)}\left[\mathrm{\Theta }_{AB}\frac{I_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}+\frac{K_0(\lambda _Fr_I)}{I_0(\lambda _Ir_I)}\right]\right)^2,$$
(53)
and for three dimensions,

$$\tau =\left(\frac{i_0(\lambda _Fr_0)}{\mathrm{\Theta }_{AB}i_0(\lambda _Fr_0)+k_0(\lambda _Fr_0)}\left[\mathrm{\Theta }_{AB}\frac{i_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}+\frac{k_0(\lambda _Fr_I)}{i_0(\lambda _Ir_I)}\right]\right)^2,$$
(54)
since $`\lambda _F>\lambda _I`$ we can verify that $`T_1(0)>T_2(0)`$, this means $`\tau >1`$. So, the temperatures ratio between the region with intruder and the region without intruder are different. In our model this lack of equipartition is due to a difference between the collisional dissipation rate related to the particle sizes. In the elastic limit $`e1`$ the energy equipartition is restored $`\tau 1`$. In Fig. 2, we present the qualitative behavior of $`\tau `$ with the size ratio $`\varphi =r_I/r_F`$, for different values of the coefficient $`e`$. The granular temperature difference increases with $`\varphi `$ and depends on $`e`$. We can see that $`\tau `$ is nearly constant and very close to unity when $`e=0.99`$.

Recently, this quantity was directly measured in experiments performed by Feitosa and Menon. They observed that energy equipartition does not generally hold for a binary vibrated granular system. They reported that the ratio of granular temperatures depends on the ratio of particle mass densities. Also in fluidized binary granular mixtures the breakdown of the energy equipartition was observed experimentally and described theoretically in the framework of the kinetic theory

## IV Thermal expansion

The definition of the granular temperature allows us to establish a statistical mechanics description of the granular medium. The change of mean energy of the system is basically due to a mechanical interaction with their external parameters (e.g., the amplitude $`A_0`$ and the frequency $`\omega _0=2\pi f`$ of vibration, the volume of the system $`V`$, and the pressure $`p`$). Theoretically we can derive the energy relaxation to the steady state for a driven granular medium. The system increases its energy as a result of external driving while its decreases its energy by dissipation. The work $`W`$ done to change the volume of the system from $`V`$ to a certain quantity $`V+dV`$ is equal to the change of its mean energy and its related to the mean pressure and volume by $`dW=pdV+Vdp`$. From the definition of granular temperature, the change of the granular temperature depends on the mean kinetic energy of the particles. A volume change $`dV`$ is related to a temperature change $`dT_g`$ by the equation of state (6).

We can express $`V`$ as a function of $`T_g`$ and $`p`$, $`V=V(T_g,p)`$. Thus given infinitesimal changes in $`T_g`$ and $`p`$, we can write

$$\begin{array}{cc}dV\hfill & =\left(\frac{V}{T_g}\right)_pdT_g+\left(\frac{V}{p}\right)_{T_g}dp,\hfill \\ & =\alpha VdT_gk_pVdp,\hfill \end{array}$$
(55)
where $`\alpha `$ is the thermal expansion coefficient defined as

$$\alpha \frac{1}{V}\left(\frac{V}{T_g}\right)_p=\frac{1}{n}\left(\frac{n}{T_g}\right)_{p,N},$$
(56)
and $`k_p`$ is the “isothermal compressibility” defined as

$$k_p\frac{1}{V}\left(\frac{V}{p}\right)_{T_g}=\frac{1}{n}\left(\frac{n}{p}\right)_{T_g}.$$
(57)
If in a first approximation we neglect the variations of the coefficients $`\alpha `$ and $`k_p`$, we can integrate Eq.(55) and find

$$\begin{array}{cc}V(T_g,p)\hfill & =V_0\mathrm{exp}\left[\alpha \mathrm{\Delta }T_gk_p\mathrm{\Delta }p\right],\hfill \\ & V_0\left[1+\alpha \mathrm{\Delta }T_gk_p\mathrm{\Delta }p\right].\hfill \end{array}$$
(58)
From the temperature difference, Eqs.(41) and (52), we can conclude that we also have a change in the effective density ($`\rho =nm_I`$) of the granular fluid between case 1 and 2, through the thermal expansion of volume produced by the difference of temperatures, Eq.(58), thus

$$\stackrel{~}{\rho }=\rho (1\alpha \mathrm{\Delta }T_g).$$
(59)
The thermal expansion coefficient can be derived from the equation of state (6) and definition (56). The general form of the coefficient of thermal expansion is

$$\alpha =\frac{1}{T_0}C(\nu ),$$
(60)
where $`C(\nu )`$ is a correction due to the density of the system. In the dilute limit $`\nu 0`$ and $`C(\nu )1`$, and the above expression tends to the expected value for a classical gas $`\alpha =1/T_0`$. The explicit form of $`C(\nu )`$ is given in Appendix B.

## V Segregation forces

Now we propose that this density difference leads to an effective buoyancy force $`𝐟_b`$, similar to the Archimedean force

$$𝐟_b=\mathrm{\Delta }\rho V_I𝐠,$$
(61)
where $`\mathrm{\Delta }\rho =\alpha \rho \mathrm{\Delta }T_g`$, $`V_I=\frac{\mathrm{\Omega }_D}{D}r_I^D`$ is the $`D`$–dimensional volume of the intruder and $`𝐠`$ is the gravity field. The intruder also experiences a viscous drag of the granular fluid. The drag force $`𝐟_d`$ is considered to be linear in the velocity of segregation $`𝐮(t)`$, and is like the Stokes’ drag force

$$𝐟_d=6\pi \mu r_I𝐮(t),$$
(62)
where $`\mu `$ is the coefficient of viscosity of the granular fluid The state–dependent viscosity posesses the general form

$$\mu =\mu _0\sqrt{T_g},$$
(63)
where the prefactor $`\mu _0`$ is a function of the fluid particle properties, and can be calculated using a Chapman–Enskog procedure for the solution of Enskog transport equation. The explicit expressions of these prefactors are given in Appendix A.

Equations (61) and (62) express the acting forces in the segregation process

$$𝐟_{seg}=𝐟_b+𝐟_d.$$
(64)
Therefore, the equation of motion that governs the segregation process is

$$\frac{\mathrm{\Omega }_D}{D}r_I^D\rho \frac{d𝐮(t)}{dt}=\frac{\mathrm{\Omega }_D}{D}r_I^D\alpha \rho \mathrm{\Delta }T_g𝐠6\pi \mu r_I𝐮(t).$$
(65)
Now we suppose the granular system contained between two large parallel plates perpendicular to the gravitational field. We take the reference frame positive in the upward vertical direction. Arranging terms in Eq.(65) we find the following differential equation

$$\frac{du(t)}{dt}=\alpha \mathrm{\Delta }T_gg\frac{6\pi D\mu (\varphi r_F)^{1D}}{\mathrm{\Omega }_D\rho }u(t),$$
(66)
where we have expressed the intruder’s radius as function of the size ratio dependence $`r_I=\varphi r_F`$, and the solution of this differential equation is the rise velocity of the intruder

$$u(t)=\frac{\alpha \mathrm{\Delta }T_ggt_0}{\varphi ^{1D}}\left[1\mathrm{exp}\left(\varphi ^{1D}\frac{t}{t_0}\right)\right],$$
(67)
where the time–scale $`t_0`$ is given by

$$t_0\frac{\mathrm{\Omega }_D\rho }{6\pi D\mu r_F^{1D}}.$$
(68)
The force balance between the drag force $`𝐟_d`$ and the buoyant force $`𝐟_b`$ gives the settling velocity $`u_s`$

$$u_s=\frac{\alpha \mathrm{\Delta }T_gt_0}{\varphi ^{1D}}.$$
(69)
The time dependent intruder height $`z(t)`$ is

$$z(t)=\frac{\alpha \mathrm{\Delta }T_ggt_0}{\varphi ^{1D}}\left[tt_0\varphi \left(1\mathrm{exp}\left(\varphi ^{1D}\frac{t}{t_0}\right)\right)\right].$$
(70)
On a qualitative level our model satisfactorily reproduces the observed phenomenology: a large intruder migrates to the top of a vibrated bed, and the rise velocity increases with the intruder size. The solutions (67) and (70) are plotted in Figs. 3 and 4. Our results resembles the experimental intruder height time evolution described in Refs. and . However the model can not describe the intermittent ascent of the intruder since we calculate the mean velocities. Using the following model parameters: mass particle density of 2.7 gcm<sup>-3</sup> (Aluminum), $`r_F=0.1`$ cm, $`e=0.83`$, $`\nu =0.34`$, $`N=5\times 10^3`$, $`g=100`$ cms<sup>-2</sup>, $`r_0=2r_I`$ and $`T_0=30\times 10^6`$ gcm<sup>2</sup>s<sup>-2</sup> we obtain that the order of magnitude of $`z(t)`$ (Fig. 4) coincides with the values reported by Cooke et al. (See Fig. 3 Ref.).

From the settling velocity $`u_s`$ (69) we show explicitly the dependence on size. This is proportional to the size ratio $`\varphi `$ and the granular temperature differences $`\mathrm{\Delta }T_g`$ which also depends on the size ratio. Its agrees with the experimental fact that the larger the radius of the intruder, the faster is the ascent, reported by Duran et al.. The plotted solution (69), Fig. 5, fits well with the experimental values from Ref. for $`\varphi >4`$, shown in Fig. 5. In this experiment, Duran et al., claim the experimental evidence of a segregation size threshold at $`\varphi _c=3.3`$, below which the intruder does not exhibit any upwards motion. Our model’s continuous aspect doesn’t allow for the existence of this threshold. We argue that this discrepancy comes from the fact that experimental measures in this regime should be very difficult to carring out. They reported to not have observed any upward motion after one hour. If we see the height profile shown in Fig. 4 for $`\varphi =2`$, we note that it is very close to zero after one minute and the slope is very low. Therefore this threshold should not exist.

## VI Conclusions

We derived a phenomenological continuum description for particle size segregation in granular media. We propose a buoyancy–driven segregation mechanism caused by the dissipative nature of the collisions between grains. The collisional dissipation rate naturally leads to a local temperature difference among the region around the intruder and the medium without intruder. In this model we proposed that the intruder’s presence develops a temperature gradient in the system which gives origin to a difference of densities. The granular temperature difference is due to the fact that the number of collisions on the surface increases with the size of the intruder, but the local density of dissipated energy diminishes. So, the region around the intruder its hotter than the region without intruder. From this temperature difference we can conclude that we have a change in the effective density of the granular fluid. This leads to an effective buoyancy force that is the responsible for the intruder’s upward movement.

In this work we made use of the tools of kinetic theory of gases to calculate the granular temperature. We observed a breakdown of the energy equipartition. And this is in agreement with other reported experiments and models. In a certain sense our theory unifies the different aspects observed in the size segregation phenomenon. Explicit solutions of the dependence of height and velocity are calculated. The geometrical effect of a segregation threshold is not supported by our model. The intruder size dependence appears naturally in our model.

Very recently it has been shown experimentally  and by computer simulations, that the convection phenomenon in granular fluids comes from the effect of spontaneous granular temperature gradients, due to the dissipative nature of the collisions. This temperature gradient leads to a density variations. The convection rolls are caused by buoyancy effects initiated by enhanced dissipation at the walls and the tendency of the grains at the center to rise. So, this segregation mechanism could be described in the hydrodynamic framework proposed in this work subject to the appropriate boundary conditions.

###### Acknowledgements.

We thank M. Alam for helpful comments on the manuscript. One of the authors (L.T.) would like to thank A.R. Lima for friendly support.

## A Transport coefficients

In this appendix the prefactors appearing in Eqs.(9) and (63) are derived. Using a Chapman–Enskog procedure for the solution of the Enskog transport equation, the transport coefficients for nearly elastic particles have been derived in Refs. and .

In $`2D`$ the thermal conductivity $`\kappa `$ is

$$\kappa =3nr_F\left(\frac{\pi }{m_F}\right)^{1/2}\left[1+\frac{1}{3}\frac{1}{G}+\frac{3}{4}\left(1+\frac{16}{9\pi }\right)G\right]T_g^{1/2},$$
(A1)
where $`G`$ is $`\nu g_0`$, $`g_0`$ is the $`2D`$ pair correlation function given in Eq.(7), and $`\nu `$ is the area fraction $`\nu =n\pi r_F^2`$. It is convenient to express Eq.(A1) introducing the prefactor $`\kappa _0`$ deffined as

$$\kappa _03nr_F\left(\frac{\pi }{m_F}\right)^{1/2}\left[1+\frac{1}{3}\frac{1}{G}+\frac{3}{4}\left(1+\frac{16}{9\pi }\right)G\right].$$
(A2)
The result (A1) takes the form

$$\kappa =\kappa _0\sqrt{T_g}.$$
(A3)
In $`3D`$ the thermal conductivity is

$$\kappa =\frac{15}{8}nr_F\left(\frac{\pi }{m_F}\right)^{1/2}\left[1+\frac{5}{24}\frac{1}{G}+\frac{6}{5}\left(1+\frac{32}{9\pi }\right)G\right]T_g^{1/2},$$
(A4)
where $`G`$ is $`\nu g_0`$, $`g_0`$ is the $`3D`$ pair correlation function given in Eq.(8), and $`\nu `$ is in this case the volume fraction $`\nu =4\pi nr_F^3/3`$. In $`3D`$ the prefactor $`\kappa _0`$ is defined as

$$\kappa _0\frac{15}{8}nr_F\left(\frac{\pi }{m_F}\right)^{1/2}\left[1+\frac{5}{24}\frac{1}{G}+\frac{6}{5}\left(1+\frac{32}{9\pi }\right)G\right].$$
(A5)
The shear viscosity $`\mu `$ in $`2D`$ is

$$\mu =\frac{1}{4}nr_F(\pi m_F)^{1/2}\left[2+\frac{1}{G}+\left(1+\frac{8}{\pi }\right)G\right]T_g^{1/2}.$$
(A6)
It is convenient to express Eq.(A6) introducing the prefactor $`\mu _0`$ defined as

$$\mu _0=\frac{1}{4}nr_F(\pi m_F)^{1/2}\left[2+\frac{1}{G}+\left(1+\frac{8}{\pi }\right)G\right].$$
(A7)
So, the result (A6) takes the form

$$\mu =\mu _0\sqrt{T_g}.$$
(A8)
In $`3D`$ the shear viscosity is

$$\mu =\frac{1}{3}nr_F(\pi m_F)^{1/2}\left[1+\frac{5}{16}\frac{1}{G}+\frac{4}{5}\left(1+\frac{12}{\pi }\right)G\right]T_g^{1/2},$$
(A9)
and the prefactor $`\mu _0`$ in $`3D`$ is defined as

$$\mu _0=\frac{1}{3}nr_F(\pi m_F)^{1/2}\left[1+\frac{5}{16}\frac{1}{G}+\frac{4}{5}\left(1+\frac{12}{\pi }\right)G\right].$$
(A10)
## B Thermal expansion coefficient

We can consider the volume of the system as a function of the granular temperature and the pressure $`V=V(T_g,p)`$. A change in the granular temperature $`dT_g`$ and the pressure $`dp`$, leads to the corresponding change in the volume $`dV`$

$$dV=\left(\frac{V}{T_g}\right)_pdT_g+\left(\frac{V}{p}\right)_{T_g}dp.$$
(B1)
As we have supposed that the pressuer of the system is more or less constant, we can approximate $`dp0`$. The increment of volume $`dV`$ with an increment of the granular temperature $`dT_g`$ is

$$dV=\left(\frac{V}{T_g}\right)_pdT_g.$$
(B2)
Thus,

$$\frac{dV}{dT_g}=\left(\frac{V}{T_g}\right)_p,$$
(B3)
or

$$\left(\frac{V}{T_g}\right)_p=\left[\left(\frac{T_g}{V}\right)_p\right]^1,$$
(B4)
and in terms of the number density $`n`$, we have

$$\left(\frac{n}{T_g}\right)_{p,N}=\left[\left(\frac{T_g}{n}\right)_{p,N}\right]^1.$$
(B5)
From the definition of the coefficient of thermal expansion Eq.(56), and from the above statement, we find

$$\alpha =\frac{1}{n}\left(\frac{n}{T_g}\right)_{p,N}=\left[\left(\frac{T_g}{n}\right)_{p,N}\right]^1.$$
(B6)
The partial derivative $`(T_g/n)_{p,N}`$ can be calculated from the equation of state (B7). In $`2D`$ the equation of state is

$$p=\frac{1+e}{2}nT_g\left[1+2\nu \frac{\left(1\frac{7}{16}\nu \right)}{\left(1\nu \right)^2}\right],$$
(B7)
where $`\nu =n\pi r_F^2`$. So, an elementary calculation leads to

$$\left(\frac{T_g}{n}\right)_{p,N}=\frac{2}{(1+e)}\frac{p}{n^2}\frac{8(\nu ^33\nu ^28\nu 8)(\nu 1)}{(\nu ^2+8)^2}.$$
(B8)
From Eq.(B6) one obtains:

$$\alpha =\frac{2}{1+e}\frac{n}{p}\frac{(\nu ^2+8)^2}{8(\nu ^33\nu ^28\nu 8)(\nu 1)}.$$
(B9)
Using the equation of state (B7) we can express $`\alpha `$ in function of the granular temperature

$$\alpha =\frac{1}{T_g}\frac{(\nu ^2+8)^2}{(\nu ^33\nu ^28\nu 8)(\nu 1)},$$
(B10)
this is

$$\alpha =\frac{1}{T_g}C(\nu ),$$
(B11)
where the correction coefficient due to the density of the system is defined as

$$C(\nu )\frac{(\nu ^2+8)^2}{(\nu ^33\nu ^28\nu 8)(\nu 1)}.$$
(B12)
For three dimensions the equation of state is

$$p=\frac{1+e}{2}nT_g\left[1+4\nu \frac{(2\nu )}{2(1\nu ^3)}\right].$$
(B13)
In a similar way we find for $`3D`$ that the coefficient of thermal expansion is

$$\alpha =\frac{1}{T_g}\frac{(\nu ^3\nu ^2\nu 1)(\nu 1)}{(\nu ^44\nu ^3+4\nu ^2+4\nu +1)},$$
(B14)
and the correction coefficient $`C(\nu )`$ in $`3D`$ is defined as

$$C(\nu )\frac{(\nu ^3\nu ^2\nu 1)(\nu 1)}{(\nu ^44\nu ^3+4\nu ^2+4\nu +1)},$$
(B15)
where $`\nu =4n\pi r_F^3/3`$.


text:
# The “action” variable is not an invariant for the uniqueness in the inverse scattering problem

## 1. Jacobi matrices of Szegö class. Direct scattering — Bernstein–Szegö type theorem

As it is well known in the theory of completely integrable systems the absolute values of a reflection coefficient have played the role of the “action” variables and the arguments of this function have meaning of the “angle” variables (in this case we think on the Toda lattice as on an integrable system) . Combining results of our previous works  (see also ) and , we give a wide set of examples where two reflection coefficients, having the same absolute values, possess completely different properties: the first one uniquely determines a Jacobi matrix of Szeqö class and the second one does not. Note that the proof of the main theorem in  (and, therefore, our result) essentially uses the analysis of  on regularization of so called Arov–singular matrix functions.

Let $`J`$ be a Jacobi matrix defining a bounded self–adjoint operator on $`l^2()`$:

(1)
$$Je_n=p_ne_{n1}+q_ne_n+p_{n+1}e_{n+1},n,$$

where $`\{e_n\}`$ is the standard basis in $`l^2()`$, $`p_n>0`$. The resolvent matrix–function is defined by the relation

(2)
$$R(z)=R(z,J)=^{}(Jz)^1,$$

where $`:^2l^2()`$ is such that

$$\left[\begin{array}{c}c_1\\ c_0\end{array}\right]=e_1c_1+e_0c_0.$$

This matrix–function possesses an integral representation

(3)
$$R(z)=\frac{d\sigma }{xz}$$

with a $`2\times 2`$ matrix–measure having a compact support on $``$. $`J`$ is unitary equivalent to the multiplication operator by an independent variable on

$$L_{d\sigma }^2=\{f=\left[\begin{array}{c}f_1(x)\\ f_0(x)\end{array}\right]:f^{}𝑑\sigma f<\mathrm{}\}.$$

The spectrum of $`J`$ is called absolutely continuous if the measure $`d\sigma `$ is absolutely continuous with respect to the Lebesgue measure on the real axis,

(4)
$$d\sigma (x)=\rho (x)dx.$$

It is natural to ask how properties of coefficients of $`J`$ are reflected on its spectral properties. One is especially interested in $`J`$’s “close” to the “free” matrix $`J_0`$ with constant coefficients, $`p_n=1,q_n=0`$ (so called Chebyshev matrix).

Let us mention that $`J_0`$ has the following functional representation, besides the general one mentioned above. The resolvent set of $`J_0`$ is the domain $`\overline{}[2,2]`$. Let $`z(\zeta ):𝔻\overline{}[2,2]`$ be a uniformization of this domain, $`z(\zeta )=1/\zeta +\zeta `$. With respect to the standard basis $`\{t^n\}_n`$ in

$$L^2=\{f(t):_𝕋|f|^2𝑑m\},$$

the matrix of the operator of multiplication by $`z(t),t𝕋`$, is the Jacobi matrix $`J_0`$, since $`z(t)t^n=t^{n1}+t^{n+1}`$.

We say that $`J`$ with absolutely continuous spectrum $`[2,2]`$ is of Szegö class if its spectral density (4) satisfies

(5)
$$\mathrm{log}det\rho (z(t))L^1.$$

###### Theorem 1.1.

Let $`J`$ be of Szegö class. Then

(6)
$$p_n1,q_n0,n\pm \mathrm{}.$$

Moreover, there exist generalized eigenvectors

(7)
$$\begin{array}{cc}& p_ne^+(n1,t)+q_ne^+(n,t)+p_{n+1}e^+(n+1,t)=z(t)e^+(n,t)\hfill \\ & p_ne^{}(n,t)+q_ne^{}(n1,t)+p_{n+1}e^{}(n2,t)=z(t)e^{}(n1,t)\hfill \end{array}$$

such that the following asymptotics hold true

(8)
$$\begin{array}{cc}\hfill s(t)e^\pm (n,t)=& s(t)t^n+o(1),n+\mathrm{}\hfill \\ \hfill s(t)e^\pm (n,t)=& t^n+s_{}(t)t^{n1}+o(1),n\mathrm{}\hfill \end{array}$$

in $`L^2`$.

Thus the eigenvectors of $`J`$ behave asymptotically as the eigenvectors of $`J_0`$ (later we make more precise statement).

The matrix formed by the coefficients of (8)

(9)
$$S(t)=\left[\begin{array}{cc}s_{}& s\\ s& s_+\end{array}\right](t),t𝕋,$$

is called the scattering matrix of $`J`$. It is unitary–valued, possesses the symmetry property $`S^{}(\overline{t})=S(t)`$ and the following analytic property: $`s(t)`$ is an outer function.

In what follows every matrix–function of the form (9) with the above listed properties is called a scattering matrix. Of course we have a good reason for this, since with every matrix $`S(t)`$ of this kind one can associate a Jacobi matrix $`J`$ whose scattering matrix (associated to $`J`$ according to Theorem 1.1) is the initial matrix–function S(t). However, $`S(t)`$, generally speaking, does not determine $`J`$ uniquely.

To clarify all above statements we need some notation and definitions. First of all for a given functions $`s_\pm `$ we define the metric

$$\begin{array}{cc}\hfill f_{s_\pm }^2=& \frac{1}{2}\left[\begin{array}{cc}1& \overline{s_\pm (t)}\\ s_\pm (t)& 1\end{array}\right]\left[\begin{array}{c}f(t)\\ \overline{t}f(\overline{t})\end{array}\right],\left[\begin{array}{c}f(t)\\ \overline{t}f(\overline{t})\end{array}\right]\hfill \\ \hfill =& f(t)+\overline{t}(s_\pm f)(\overline{t}),f(t),fL^2,\hfill \end{array}$$

and we denote by $`L_{dm,s_\pm }^2`$ or $`L_{s_\pm }^2`$ (for shortness) the closure of $`L^2`$ with respect to this new metric.

The following relations set a unitary map from $`L_{s_+}^2`$ to $`L_s_{}^2`$:

(10)
$$\begin{array}{cc}\hfill \left[\begin{array}{c}sf^+\\ sf^{}\end{array}\right](t)=& \left[\begin{array}{cc}s& 0\\ s_+& 1\end{array}\right](t)\left[\begin{array}{c}f^+(t)\\ \overline{t}f^+(\overline{t})\end{array}\right]\hfill \\ \hfill =& \left[\begin{array}{cc}1& s_{}\\ 0& s\end{array}\right](t)\left[\begin{array}{c}\overline{t}f^{}(\overline{t})\\ f^{}(t)\end{array}\right].\hfill \end{array}$$

Moreover, in this case,

$$f^+_{s_+}^2=f^{}_s_{}^2=\frac{1}{2}\{sf^+^2+sf^{}^2\}.$$

It is worth to give a scalar variant of relations between $`f^+L_{s_+}^2`$ and $`f^{}L_s_{}^2`$:

$$s(t)f^{}(t)=\overline{t}f^\pm (\overline{t})+s_\pm (t)f^\pm (t).$$

###### Theorem 1.2.

$`J`$ is a Jacobi matrix of Szegö class with the spectrum $`E=[2,2]`$ if and only if $`J`$ possesses the scattering representation, i.e.: there exists a unique matrix–function $`S(t)`$ of the form (9) (with the listed properties) and a unique pair of Fourier transforms

(11)
$$^\pm :l^2()L_{s_\pm }^2,(^\pm Jf)(t)=z(t)(^\pm f)(t),$$

determining each other by the relations

(12)
$$s(t)(^\pm f)(t)=\overline{t}(^{}f)(\overline{t})+s_{}(t)(^{}f)(t),$$

and having the following analytic properties

(13)
$$s^\pm (l^2(_\pm ))H^2,$$

and asymptotic properties

(14)
$$e^\pm (n,t)=t^n+o(1)\text{in}L_{s_\pm }^2,n+\mathrm{},$$

where

$$e^+(n,t)=(^+e_n)(t),e^{}(n,t)=(^{}e_{n1})(t),$$

with $`\{e_n\}`$ being the standard basis in $`l^2()`$.

We point out that asymptotic relations (8) and (14) are equivalent, moreover (14) directly shows that the eigenvectors of $`J_0`$ asymptotically are the eigenvectors of $`J`$.

## 2. Uniqueness and Completeness

Before we proceed with the uniqueness theorem, we show how to construct at least one $`J`$ with the given scattering matrix $`S(t)`$.

Consider the space

$$H_{s_+}^2=\text{clos}_{L_{s_+}^2}H^2,$$

and introduce the Hankel operator $`_{s_+}:H^2H^2`$,

$$_{s_+}f=P_+\overline{t}(s_+f)(\overline{t}),fH^2,$$

where $`P_+`$ is the Riesz projection from $`L^2`$ onto $`H^2`$. This operator determines the metric in $`H_{s_+}^2`$:

$$\begin{array}{cc}\hfill f_{s_+}^2=& f(t)+\overline{t}(s_+f)(\overline{t}),f(t)\hfill \\ \hfill =& (I+_{s_+})f,f,fH^2.\hfill \end{array}$$

###### Theorem 2.1.

Let $`S(t)`$ be a scattering matrix, i.e., the matrix of the form (9) with listed properties. Then the space $`H_{s_+}^2`$ is a space of holomorphic functions with a reproducing kernel. Moreover, the reproducing vector $`k_{s_+}`$:

$$f,k_{s_+}=f(0),fH_{s_+}^2,$$

is of the form

(15)
$$k_{s_+}=(I+_{s_+})^{[1]}\text{1}:=\underset{ϵ0^+}{lim}(ϵ+I+_{s_+})^1\text{1}\text{in}L_{s_+}^2.$$

Put $`K_{s_+}(t)=k_{s_+}(t)/\sqrt{k_{s_+}(0)}`$. Then the system of functions $`\{t^nK_{s_+t^{2n}}(t)\}_n`$ forms an orthonormal basis in $`L_{s_+}^2`$. With respect to this basis, the multiplication operator by $`z(t)`$ is a Jacobi matrix $`J=J[s_+]`$ of Szegö class. Moreover, the initial $`S(t)`$ serves as the scattering matrix–function, associated with given $`J`$ by Theorem 1.2, and the relations

$$^+(e_n)=t^nK_{s_+t^{2n}}(t)$$

determine uniquely corresponding Fourier transforms.

Let us fix the notation $`J[s_+]`$ for the Jacobi matrix associated with $`S(t)`$ by Theorem 2.1. On the other hand, the system of functions $`\{t^nK_{s_{}t^{2n}}(t)\}_n`$ forms an orthonormal basis in $`L_s_{}^2`$, and we can define a Jacobi matrix $`\stackrel{~}{J}=J[s_{}]`$ by the relation

$$z(t)\stackrel{~}{e}^+(n,t)=\stackrel{~}{p}_n\stackrel{~}{e}^+(n1,t)+\stackrel{~}{q}_n\stackrel{~}{e}^+(n,t)+\stackrel{~}{p}_{n+1}\stackrel{~}{e}^+(n+1,t),$$

where $`\{\stackrel{~}{e}^+(n,t)\}`$ is the dual system to the system $`\{t^nK_{s_{}t^{2n}}(t)\}`$ (see (12)), i.e.:

$$s(t)\stackrel{~}{e}^+(n1,t):=\overline{t}^{n+1}K_{s_{}t^{2n}}(\overline{t})+s_{}(t)t^nK_{s_{}t^{2n}}(t).$$

None guarantees that operators $`J[s_+]`$ and $`J[s_{}]`$ are the same (see beginning of the next section). However, if $`J[s_+]=J[s_{}]`$, then the uniqueness theorem takes place.

###### Theorem 2.2.

A scattering matrix $`S(t)`$ determines a Jacobi matrix $`J`$ of Szegö class in a unique way if and only if the following relations take place

(16)
$$s(0)K_{s_\pm }(0)K_{s_{}t^2}(0)=1.$$

Of course it is hard to check identities, especially using computer simulation, but, in fact, (16) has a specific approximating meaning.

Let us return to the matrix $`J[s_{}]`$. According to (13) the space $`^+(l^2(_+))`$ is a subspace of $`L_{s_+}^2`$ consisting of holomorphic in $`𝔻`$ functions. In the given case we can even specify this space in the form

(17)
$$^+(l^2(_+))=\widehat{H}_{s_+}^2:=\{fL_{s_+}^2:sfH^2\}.$$

Every function $`f`$ from $`H_{s_+}^2`$ possesses the property $`sfH^2`$, but this means only inclusion:

(18)
$$\widehat{H}_{s_+}^2H_{s_+}^2$$

The meaning of (16) is that every function from $`\widehat{H}_{s_+}^2`$ can be approximated by functions from $`H^2`$ in $`L_{s_+}^2`$–norm ($`H^2`$ is dense in $`\widehat{H}_{s_+}^2`$). In fact, it is enough to prove that we can approximate just two functions $`\stackrel{~}{e}^+(0,t)`$ and $`\stackrel{~}{e}^+(1,t)`$ from $`\widehat{H}_{s_+}^2`$. This would guarantee (16), completeness and the uniqueness theorem.

## 3. The result

Our first remark is almost evident. Let us pick an analytic scattering matrix (all entries are in $`H^{\mathrm{}}`$). For example,

$$S(t)=\left[\begin{array}{cc}\frac{1+\mathrm{\Delta }}{2}& \frac{1\mathrm{\Delta }}{2}\\ \frac{1\mathrm{\Delta }}{2}& \frac{1+\mathrm{\Delta }}{2}\end{array}\right](t),$$

where $`\mathrm{\Delta }(t)=\overline{\mathrm{\Delta }(\overline{t})}`$ is a symmetric inner function. Then the analyticity of $`s_+`$ implies that $`_{s_+}=0`$ , $`H_{s_+}^2=H^2`$ and $`t^nK_{s_+t^{2n}}(t)=t^n`$, $`n0`$. Thus,

$$J[s_+]|l^2(_+(1))=J_0|l^2(_+(1)),$$

where $`_+(m)=\{n:nm\}`$. The analyticity of $`s_{}`$, in its turn, implies that

$$J[s_{}]|l^2(_{}(1))=J_0|l^2(_{}(1)),$$

$`_{}(m)=\{n:nm1\}`$.

Thus, the uniqueness would imply that $`J`$ associated to given $`S(t)`$ coincides with $`J_0`$ up to three coefficients: $`p_0,q_0,q_1`$. For fun we give an exact formula for $`S(t)`$ related to this case:

(19)
$$S(t)=\frac{1}{\varphi (t)}\left[\begin{array}{cc}\psi _{}(t)& p_0(1t^2)\\ p_0(1t^2)& \psi (t)\end{array}\right],$$

where $`\varphi (t)=(1q_0t)(1q_1t)p_0^2t^2`$, $`\psi (t)=(1q_0t)(q_1t)+p_0^2t^2`$ and $`\psi _{}(t)=t^2\overline{\psi (t)}`$. Note that this matrix has no bound states ($`s(t)`$ is holomorphic in $`𝔻`$) if

$$\left[\begin{array}{cc}q_1& p_0\\ p_0& q_0\end{array}\right]1.$$

Therefore, once the scattering matrix is not a rational function of degree two, then, for sure, the holomorphic $`S(t)`$ does not determine $`J`$ uniquely.

The main result of this note is as follows.

###### Theorem 3.1.

Given $`S(t)`$ from $`H^{\mathrm{}}`$, one can find an inner function $`\mathrm{\Phi }(t)=\overline{\mathrm{\Phi }(\overline{t})}`$ such that

$$S_\mathrm{\Phi }(t)=\left[\begin{array}{cc}s_{}\mathrm{\Phi }& s\\ s& s_+\overline{\mathrm{\Phi }}\end{array}\right](t),t𝕋,$$

determines a Jacobi matrix of Szegö class uniquely.

Let us point out that we did not change the transmission coefficient $`s`$. Moreover the new reflection coefficient $`s_{}\mathrm{\Phi }`$ is still analytic one, so, the corresponding Jacobi matrix $`J`$ coincides with $`J_0`$, when both are restricted on $`l^2(_{}(1))`$.

Theorem 3.1 is the direct consequence of the main result of  and the following two lemmas.

###### Lemma 3.2.

Let the reflection coefficient $`s_{}H^{\mathrm{}}`$, and $`s_{}(\zeta )=a\zeta +`$, i.e.: $`s_{}(0)=0`$. Then (16) holds if

(20)
$$\left[\begin{array}{cc}1& s_{}\\ 0& s\end{array}\right]\left[\begin{array}{cc}1& \overline{t}\\ \overline{t}& 1\end{array}\right]=\underset{n\mathrm{}}{lim}\left[\begin{array}{cc}s& 0\\ s_+& 1\end{array}\right]\left[\begin{array}{cc}f_n& \overline{t}g_n\\ \overline{t}f_n(\overline{t})& g_n(\overline{t})\end{array}\right]$$

with a suitable choice of sequences of functions $`f_n,g_n`$ from $`H^2`$, that is, $`\frac{1+s_{}/t}{s}H_{s_+}^2`$ and $`\frac{1+ts_{}}{s}H_{s_+t^2}^2`$.

###### Proof.

Since $`s_{}`$ is holomorphic, we have $`_s_{}=0`$. Moreover, even $`_{s_{}t^2}\text{1}=a\text{1}`$, since $`s_{}(0)=0`$. Therefore,

(21)
$$K_s_{}=\text{1}\text{and}K_{s_{}t^2}=\frac{\text{1}}{\sqrt{1+a}}.$$

We work with the second vector. Because of

(22)
$$\begin{array}{cc}\hfill \left[\begin{array}{cc}s& 0\\ s_+& 1\end{array}\right]\left[\begin{array}{c}f^+(t)\\ \overline{t}f^+(\overline{t})\end{array}\right],\left[\begin{array}{cc}1& s_{}\\ 0& s\end{array}\right]\left[\begin{array}{c}1\\ \overline{t}\end{array}\right]=& \left[\begin{array}{cc}s& 0\\ 0& \overline{s}\end{array}\right]\left[\begin{array}{c}f^+(t)\\ \overline{t}f^+(\overline{t})\end{array}\right],\left[\begin{array}{c}1\\ \overline{t}\end{array}\right]\hfill \\ \hfill =& 2(sf^+)(0),f^+H^2,\hfill \end{array}$$

and (20), the vector $`\frac{1+s_{}/t}{s\sqrt{1+a}}`$ belongs to $`H_{s_+}^2`$ and it is collinear with the reproducing kernel of this space. Since its norm is one, we get

$$K_{s_+}(\zeta )=\frac{1+s_{}(\zeta )/\zeta }{s(\zeta )\sqrt{1+a}}.$$

Putting here $`\zeta =0`$, we get $`s(0)K_{s_+}(0)K_{s_{}t^2}(0)=1`$. The proof of the second identity is very similar. ∎

###### Lemma 3.3.

Let the following approximation hold true

(23)
$$\left[\begin{array}{c}\overline{t}s_{}\\ \frac{ss(0)}{t}\end{array}\right]=\underset{n\mathrm{}}{lim}\left[\begin{array}{c}su_n\\ P_+s_+u_n\end{array}\right],u_nH^2$$

(we still assume that $`s_{}(0)=0`$). Then (20) is also true.

###### Proof.

We put

(24)
$$\overline{t}f_n(\overline{t})=s(0)\overline{t}+\overline{t}u_n(\overline{t})P_{}s_+u_n$$

and $`g_n=f_n(1t^2)u_n`$. In this case

$$\left[\begin{array}{cc}f_n& \overline{t}g_n\\ \overline{t}f_n(\overline{t})& g_n(\overline{t})\end{array}\right]=\left[\begin{array}{cc}f_nu_n& tu_n\\ \overline{t}u_n(\overline{t})& f(\overline{t})u(\overline{t})\end{array}\right]\left[\begin{array}{cc}1& \overline{t}\\ \overline{t}& 1\end{array}\right].$$

In fact, we are going to show that

(25)
$$\left[\begin{array}{cc}1& s_{}\overline{t}\\ 0& s\overline{t}\end{array}\right]=\underset{n\mathrm{}}{lim}\left[\begin{array}{cc}s& 0\\ s_+& 1\end{array}\right]\left[\begin{array}{cc}f_nu_n& u_n\\ \overline{t}u_n(\overline{t})& \overline{t}\{f_n(\overline{t})u_n(\overline{t})\}\end{array}\right].$$

Let us note that the equality in the second column of (25) is just a direct consequence of (23) and (24). Therefore, we need to check the equality in the first column. In the first entry of this column we have

$$\begin{array}{cc}\hfill \underset{n\mathrm{}}{lim}s(f_nu_n)=& \underset{n\mathrm{}}{lim}\{sf_ns_{}\overline{t}\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{ss(0)+su_nsP_+\overline{t}s_+(\overline{t})u_n(\overline{t})s_{}\overline{t}\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{ss(0)+s_{}\overline{t}s\overline{s_+}\overline{t}u_n(\overline{t})+sP_{}\overline{t}s_+(\overline{t})u_n(\overline{t})s_{}\overline{t}\}.\hfill \end{array}$$

Using the unitary property of $`S`$–matrix and (23), we get

$$\begin{array}{cc}\hfill \underset{n\mathrm{}}{lim}s(f_nu_n)=& \underset{n\mathrm{}}{lim}\{ss(0)+s_{}\overline{s}\overline{t}u_n(\overline{t})+sP_{}\overline{t}s_+(\overline{t})u_n(\overline{t})\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{ss(0)+s_{}\overline{s_{}}+s(\overline{s}s(0))\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{|s_{}|^2+|s|^2\}=1.\hfill \end{array}$$

For the second entry we have

$$\begin{array}{cc}\hfill \underset{n\mathrm{}}{lim}\{s_+(f_nu_n)+\overline{t}u_n(\overline{t})\}=& \underset{n\mathrm{}}{lim}\{s_+s(0)s_+P_+\overline{s_+}\overline{t}u_n(\overline{t})+\overline{t}u_n(\overline{t})\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{s_+s(0)+s_+P_{}\overline{s_+}\overline{t}u_n(\overline{t})+(1|s_+|^2)\overline{t}u_n(\overline{t})\}\hfill \\ \hfill =& \underset{n\mathrm{}}{lim}\{s_+s(0)+s_+(\overline{s}s(0))+s\overline{s}\overline{t}u_n(\overline{t})\}\hfill \\ \hfill =& s_+s(0)+s_+(\overline{s}s(0))+s\overline{s_{}}=0.\hfill \end{array}$$