Christelle Garrouste
email: christelle.garrouste@univ-orleans.fr

Cyrille Piatecki

Yvan Stroppa

Matching dynamics and optima in a multi-agents labor market setting

One of the greatest difficulties meet on the labor market is that it ensures that all the people who want to work can find a job. But it is not the only one. A second difficulty consists to find a social arrangement which approaches as much as possible the social optimum without forgetting that the optimum can be reached with the rejection of certain candidacies.

It seems impossible to solve a social optimum problem of the size of the french labor market because this one is furthermore constituted of 26 million job-seekers and that no computer is capable of managing the matching problem which, surprisingly is not of an excessive complexity -it is of order n 3 -, but which quickly collides with the computational capacity of machines, however powerful they may be.

It's the reason why it is essential to study restricted size markets, that is markets of less than 10000 job-seekers, to observe how practical solutions which can be organized as, for example, a geographical progressive balkanization either geographic or by skill level, leads to a departure from the Pareto optimum.

In this context, this paper suggests to compare the level of efficiency of a decentralized actually practicable solution with the optimum which would be if only it was possible to collect all necessary information to implement and compute it.

Introduction

Labor markets are notably inefficient. There is always a huge discrepancy between a more or less important reservoir for tenders unfilled jobs and a significant unemployment. For instance, for France the DARES publishes a quarterly indicator of the tension of the labor market -see the figure 1 -which demonstrates that since 1998, and certainly earlier, the 74 professional families 1 have'nt displayed any equilibrium between supply and offer. As a high level of unemployment has always been considered as dangerous for the social status quo, after one century of intervention it is difficult to identify the short term reasons of the discrepancy. If in the long term, one knows that the two main sources come from the weakness of the activity and the impossibility to anticipate which type of formation will be asked by firms, in the short term, it is less clear what are the main sources of the inefficiency.

The reason is elementary: at least in Europa, there is no more countries purely managed through the market or through a socialist organization and even if one can see a true tendency to let the market regain the upper hand of the plan, there are good reasons to keep a status quo, even if one can not say right now what will be, in the end, the role between centralised allocation of jobs and the market. Now, since in 1908, Barone -see [START_REF] Barone | The ministry of production in the collectivist state[END_REF] -has convincingly argued that a collectivist managed economy achieves, as the market, an efficient allocation of resources and [START_REF] Arrow | Existence of an equilibrium for a competitive economy[END_REF], through the two famous welfare theorems, have shown under which conditions the market is, the Pareto criterion -voir [START_REF] Pareto | Manual of Political Economy[END_REF] -is the core description of what would be an intelligent centralized economy.

Unfortunately, even if, from the theoretical point of view, it is essential, one can doubt, without reopening the Hayekian debate -see von [START_REF] Hayek | The Road to Serfdom[END_REF] -, that a planer managing an organization constituted of agents made of flesh and bones could achieve even an ε-Pareto optimality2 .

Indeed, there is a lot of reasons not to achieve or simply to approach a Pareto optimality. One can begin by the [START_REF] Niskanen | Bureaucracy and Representative Government[END_REF] objection that any administration plays for itself -see also Buchanan et al. [1980] -, not necessarily for the collective well-being. So, it can renounce to strictly find a Pareto-optimum for an other objective less dedicated to the collective well-being.

Secondly, there is the million of equations debate, which run around the fact that since in a socialist economy there is not a self-interest mechanism there is no way to drive the exchange process to an equilibrium. The tentative answer was given as soon as [START_REF] Taylor | The guidance of production in the socialist state[END_REF] that one must not mistaken the property of the goods with the process since, coming back to Barone, the socialist economy and the capitalist economy are formally described by the same theoretical system. Two main themes were raised : the first one was about the nature of prices when Lange [1936Lange [ -1937] ] made the distinction between prices as exchange ratios and prices as terms in which alternatives are offered. This distinction open the road to the linear programming approach to the optimum, first proposed in the soviet era by Kantorovich in 1939 but which was accessible in occident only after its translation in [START_REF] Kantorovich | The mathematical method of production and planning[END_REF]. That was to say that if the central planer collects all the information, the problem is reduced to solve for prices and quantities by a simple computer program. Unfortunately, in those times, computers were realities only in the brain of a small number of scientists who do not even counted on the fingers of the hand. Therefore, no one could really discuss the time required to complete the calculations required for planning. It's only with [START_REF] Hartmanis | On the computational complexity of algorithms[END_REF] that the idea to measure the time and space required to achieve a computation emerged. Later, [START_REF] Klee | How good is the simplex algorithm? In Inequalities III (Proceedings of the Third Symposium on Inequalities held at the University of California[END_REF] have shown that, in the worst cases, the Dantzig simplex algorithm is of exponential complexity, that is to say that when the dimension of the domain is increased by one, the time necessary to achieve the optimum is lined. The consequence of this is that, in some cases, it's pure utopia to try to achieve an optimum for a huge dimension problem3 , and this is amplified if one search to set aside the hypothesis of absolute divisibility of goods which can be solved only through combinatorial algorithms whose algorithmic complexity is higher4 .

Of course, since the heroic times of the invention of linear programming in two of the most gifted mathematicians of the XXth century -Kantorovitch and Dantzig -, some more efficient interior points algorithms have been conceived which reach solutions in polynomial time.

But, if one follows [START_REF] Shalizi | In soviet union, optimization problem solves you[END_REF], even if there is now some good commercial linear programming package who could handle a 12 million variables on a desktop computer in a few minutes5 , the problem remains. This 12 million number is not randomly chosen : it appears in [START_REF] Nove | The Economics of Feasible Socialism[END_REF] as the number of identifiable distinct products in the world -it is certainly today an under-evaluation6 . But, as shown brilliantly in 1959 by Debreusee [START_REF] Debreu | Theory of Value: An Axiomatic Analysis of Economic Equilibrium[END_REF] -, argued this doesn't take in account neither the place where goods are produced, neither the place where they are consumed -to forget this simple fact, is to forget the transport services which must also be taken into account -, the time of disposal -one knows that time has a value -, the decided incorporated qualities -think simply of planed obsolescence -and the uncertainties attached to the states of Nature when it will be used -an the value of an umbrella under clear weather is not the same as the value of an umbrella under rain.

How many goods actually exist if one take into account those apparently safe little deviations from reality. Certainly a number so huge that no one could be comfortable with it and certainly no computers present or to come -even in taking in account Moore's law.

One can not hope that markets constitute a solution because they encounter the same difficulty. But, as stated by [START_REF] Axtell | The complexity of exchange[END_REF], Arrow-Debreu equilibria are sufficiently difficult to compute that the Walrasian picture of market behavior is simply not plausible7 .

Of course those are arguments that are easily amenable to static computation, but one can also speak of the necessity to collect information through processes providing strong enough incentives to warrant it's honest transmission. Then one must construct an institution efficiently organized to transform the information in a computable form.

All this take time, and more than that, gathering information is a sequential rather than a synchronic process. So, must one wait to have collected all the information to take a decision and let the system dy in waiting. In short, computing a Pareto-efficient equilibrium is accessible to an omniscient demiurge, not to a human planer. And Paretoefficiency is by now certainly the criterion by which every economic system is evaluated.

In what follows, we will study highly stylized model of a labor market. Those type of markets are characterized by a two side matching between labor demand and supply which a specific case of Matching theory, that is a mathematical account of behavior that is at the basis of many economics and computer science assignments. Stable matching problems are the most prominent models because they capture the aspect of rationality and distributed control that is inherent in most optimization problems [START_REF] Hoefer | Local matching dynamics in social networks[END_REF]). Examples of two-sided (or pairwise) matching models include the assignment of workers to jobs [START_REF] Hoefer | Local matching dynamics in social networks[END_REF], [START_REF] Arcaute | Social networks and stable matchings in the job market[END_REF] and [START_REF] Kelso | Job matchings, coalition formation and gross substitute[END_REF]), organs to patients [START_REF] Roth | Pairwise kidney exchange[END_REF]), general buyers to sellers (Lauermann, 2013) and stable marriage and stable roommate problems [START_REF] References | Autonomous actor positioning in wireless sensor and actor networks using stable-matching[END_REF], [START_REF] Goemans | Market sharing games applied to content distribution in ad-hoc networks[END_REF] and [START_REF] Mathieu | Self-stabilization in preference-based systems[END_REF]).

Initially applied to markets, the classic static Walrasian pairwise matching models [START_REF] Becker | A theory of marriage: Part i[END_REF]), assert the existence of a positive assortative matching (PAM) when productively complementary agents are assorted based upon their expected abilities. Recent developments show that the PAM still holds with random matching and search frictions [START_REF] Atakan | Assortative matching with explicit search costs[END_REF] and [START_REF] Shimer | Assortative matching and search[END_REF]).

While the literature on static matching has brought many important insights into non-cooperative games for decentralized markets, the need for the development of dynamic matching models is now widely recognized. This need is driven by the existence in the real world of recursive spells of games with evolving individual characteristics, as demonstrated by [START_REF] Anderson | Dynamic matching and evolving reputations[END_REF] 8 .

In this paper, we consider a dynamic model of stable matching in which players are rational agents seeking to be matched to their optimal partner, chosen within a social network, without any central coordination. The ordinary stable room-mates and marriage problems assert that each player knows the complete player set and can be matched arbitrarily. However, this assumption does not hold for all types of two-sided matchings. For instance, in the case of workers to jobs assignments, any worker can not be assigned arbitrarily to any job. There are restrictions in terms of knowledge and information that will hinder some players to be matched and that will enable certain players to match up faster.

To account for these matching barriers and frictions, our model builds upon the works by [START_REF] Arcaute | Social networks and stable matchings in the job market[END_REF] and [START_REF] Hoefer | Local matching dynamics in social networks[END_REF]. They developed job-market games in which they considered locally stable matchings. Their models are special cases of stable marriage, where firms strive to hire workers. The main idea is that social links exist only among workers and that each firm can match to k workers, but each worker can only match to one firm.

While we start from similar assumptions [START_REF] Arcaute | Social networks and stable matchings in the job market[END_REF] and [START_REF] Hoefer | Local matching dynamics in social networks[END_REF], we go one step ahead by generating a computerized optimization algorithm that simulates the spatial dynamics of the convergence process. This unique and innovative computer software enables us to test empirically our theoretical assumptions and to generate descriptive statics about the matching process.

This paper proceeds with a description of the model (section 1). Then we will discuss the interaction space (section 2), before the institutional arrangement preluding to the matches (section 3) and finally one will discuss the simulation results (section 4).

Description of the model

One start with two populations H and I respectively indexed by h and i -H = n and I = m.

The matching problem

Within each population, the agents are randomly assigned an individual tag t h for an H-agent and t i for an I-agent. The tag, which characterises qualitatively each agent, is randomly determined to ensure that each agent within a population is associated to a unique tag: t h ∈ {1, . . . , n} and t i ∈ {1, . . . , m}.

Then a reservation tag -t r h for the agents belonging to the H population and t r h for the agents belonging to the I population -are drown in the sets

[t h -ε h , t h + μ h ] for the H-agents [t i -ε i , t i + μ i ]
for the I-agents and according to which agents accept a match9 if t r h ≤ t i for the H-agent h to accept the pairing with the agent i t r i ≤ t h for the I-agent I to accept the pairing with the agent h

One assumes that the reservation tags are not too far of the own tag of the agents, because it seems natural to reduce the possibility of matching among agents that are too different from one another.

Respectively, if the tag of the other agent fall inside the acceptation set A of an agent the later one accepts the proposal. It is only when two agents accept a pairing that the match is possible -according to the classical voluntary exchange rule. On the contrary, if at least one tag does'nt enter in its acceptation set the match fails and both agents h, and respectively i, will have to re-initialize their search.

The demiurge arrangement 10

Before going further, one must note that, as creator of the electronic world, the authors of the simulation are in a very specific situation : they are omniscient and can calculate the best possible matches.

Whenever a h agent accepts a match he has a surplus S h ∈ [0, ε h + μ h ], and whenever an agent i accepts a match he has a surplus

S i ∈ [0, ε i + μ i ]. So, a (h/i)-match generates a surplus S hi = S h + S i . Define δ hi = 1 if h and i match 0 in the alternate case
The planer then to find the δ hi which satisfy the

O 0 = max {δ hi ,h∈H,i∈I} H h=1 I i=1 S hi δ hi under the H + I constraints for i = 1, . . . , I : H h=1 δ hi = 1 for h = 1, . . . , H : I i=1 δ hi = 1 ∀h ∈ H, ∀i ∈ I, δ hi ∈ Z
In a multi-agents simulation, as long as H +I is not too much high11 , one can compute O 0 that is the Pareto-optimal overall surplus index as the sum of all total surpluses associated with the matches. This number will be use as benchmark to evaluate the various human organization of the labor market.

As an example, the T hi table exhibited in the figure 2 gives the aggregated surpluses of all the pair (h, i) up to the restriction that the matches be accepted by the two sides. So, we can evaluate the efficiency of the spontaneous order which emerges from the two-sided search. But, one must also notice that this procedure can be such that some agents stay unmatched. For instance, the T hi table gives such a Pareto-optimal matching set which generates an overall surplus O 0 of 462, but let the agent 7 ∈ H and the agent 5 ∈ I unmatched. Let us begin by describing the operating conditions in which the simulations of the demiurge were realized. First, as regards the hardware, the software turned on a HP Z830 TM 16 core 2 Ghz which is a rather powerful and quick machine. Of course, in an idealized world, the demiurge would have been able to have access to the chinese, right now, the most powerful computer in the world -3 120 000 core. There are two types of LP-softwares : the non-commercial ones and the commercial ones. For pure doctrinal reasons, one has began by trying to resolve the some random generated situation of matching generated according to the description of the model given in the section 1, with free softwares mainly with LP-solve, because despite it's public domain licence it is a celebrated very versatile soft, running in many platforms and callable nearly from any program. The main observation in running the diverse soft is that they are time consuming. This is not surprising since the matching algorithm is of order O(n 3 ) which is harder than a less demanding O(n 2 ) algorithm, but stay in the feasible set as long as n is not too huge 12 . At least, it is not of exponential complexity, but one can easily understand that it may force to kneel down the fastest computers. It is however necessary to note that the analysis of the complexity does not reveal the problem of the size of the memory to realize the calculations. All programs, without exception, free or not, met a computational barrier as soon as n = H = I was set to 7500 and this barrier is in fact the one encountered for the commercial softwares since, for LP-solve, it was 4500. One knows that it is also possible to solve the problem using the Hungarian algorithm which is equivalent to the linear programming approach, but it the problem was of the same order. Nevertheless, it seems that CPLEX is the fastest of all the programs one has been able to use 14 .

1
All this means that, without summoning the most powerful available computers in the world, it is likely not possible to calculate the optimal matching of a labor market as the French market. And in all cases, if the the Pareto solution to the French labor market was computable, one necessarily would ask for the European labor market.

Thus it seems necessary to turn towards a second-best solution. As it will be shown in a later section, there are multiple options. The simplest one which emerged was to design an greedy algorithm with began by randomly order both sets then to propose the highest generating surplus match to the first of the H-agents and so on until no possible match remains.

In all the cases, it should be understood that in the actual world the computational dimension is one of the dimensions of the optimal matching problem and perhaps not the more striking. An employment agency has to ask the candidates to attend to a rendez-vous with one staff member. This one would have to listen to them, estimate 12 For n = 7000, the order is of 343 M operations. 13 This graph may seems weird since expect some nearly linear curves. One must realize that there is no command in those softwares that gives the exact computation time and that one must rely on the human attention to the end of the computation which are heavily soiled by errors.

14 But CPLEX fails to compute above n = 3000.

them and, in one way or an other, he will have to apply a procedure allowing him not only to attribute a tag to the the agent under inspection, but also a reservation tag. Then, a new rendez-vous with the agent turns out to be necessary to ensure that the administrative attributed tags are identical with the those that the agent appropriate to himself15 . All those operations are time consuming and the human computer, lost in the management of the necessary interaction between the administrative and the private agents, cannot manage it quickly. For all those reasons the Pareto optimal matching could not be really obtained in the actual world.

The problem with the interaction space

Actually, the agents doesn't wait in their ivory tour. Specifically, the H-agents -the laborers -move to meet a potential match, when the I-agents -the firms -wait for the first ones. That is not to say that firms have nothing to do to hire a new employee. They must make their offers public. On postulates that it is because of the Nature of the labor market that they do not move : they need not to use a head hunter because on an a priori basis there are enough laborers to fulfill all jobs. So this is an imperfect information situation for the both sides : the H-agents do not know where are the offers, the I-agents does not know where are the potential matches.

One must note that, whatever its metaphoric interpretation, the space has always been present in the search literature and more specifically, in the economic literature. Setting aside the land use models "à la" Thunen [1966] and the industrial localization models "à la" [START_REF] Hoteling | Stability in competition[END_REF], because in those two approaches to economics modeling there was no search of any information -the former was trying to explain the soil rents, and the later, the best place to locate a shop, a plant or any economic activity -, it was necessarily present in the path breaking work of [START_REF] Koopman | Search and screening[END_REF], since Koopman was mainly interested in hunting submarines. His work has been extended by [START_REF] Guenin | Optimal allocation of effort: An extension of the Koopman basic theory[END_REF] -see also [START_REF] Stone | Theory of Optimal Search[END_REF] for a recension of all the work done on the subject.

When [START_REF] Stigler | The economics of information[END_REF] imported the technology in labor economics, it was also the case because he postulated that the agents were unaware of which are the posted prices but not of where they were located, and Stigler postulated implicitly that any location could be reached from any location for a constant cost. That is all information locations are situated on the same circle centered on the searching agent and this is not acted upon if the agents choose a reservation price strategy "à la" [START_REF] Lippman | The economics of job search: A survey[END_REF] rather than a fixed sample strategy discussed by Stigler or even a mixed-one "à la" [START_REF] Morgan | Optimal search[END_REF] or [START_REF] Benhabib | Job search: The choice of intensity[END_REF].

With [START_REF] Sakoda | Minidoka: An Analysis of Changing Patterns of Social Interaction[END_REF] -see also [START_REF] Sakoda | The checkerboard model of social interaction[END_REF][START_REF] Schelling | Models of segregation[END_REF] 16 , the agents are searching for the best location possible on a chess or a checker-board according to x y z { Figure 4: Neighborhood and Small Worlds the proportion of non desired persons in their neighborhood. The space is then divided in cells and the agents can jump from a cell to an other. This create a discrete space from a continuous one.

In the same spirit, the spatial prisoners' dilemma, initially introduced by Nowak and May [1992], works through a cells organized space but, this time, with fixed agents who play only with their neighbours17 and adapt their strategy on a local information basis. This simple structure has a natural drawback: to guaranty that some agents are not in a singular situation, which set the number of their neighbours to a lower number that the most part of the other agents -a side effect -, it was necessary to dot the space with a property which make him topologically equivalent to a torus. Later with the introduction of the small world network of [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF], the space was structured as a network in which most nodes can be reached from every other by a small number of hops or steps.

In the contrary of all those specific spaces, the space here is a continuous one18 and the agents can a priori19 explore each atom of the space -i.e. each pixel. Even if, in this case, this is not truly necessary, we have imposed a basic torus structure which is no more than a topological transform of a rectangle of whichever dimension20 .

In this space, firms are thrown on fixed locations at random, as well as workers to the restriction that as soon as one simulation begins, workers move toward one objective according to the institutional arrangement which describe the market. In this paper, one will study the institutional arrangement described in the next section, under the voluntary exchange hypothesis.

The nature of the space is but one of the conditions preluding to the matches. A second important property of the space is the number of inhabitant it may welcome.

Figure 5: A Torus Along all simulation, the space has been set to 800 × 500 = 4000 pixels. As an electronic agent occupy 50 pixels, one can have at most 8000 electronic agents of one or the other type in the world. In fact, as each firm offers only one job, workers can spatially recover the firms and one reach a maximum number of 16 000 agents. But, let us keep in mind the number of 8000, because this number corresponds to a saturated world where any movement is impossible. This leads to the agents' density described in figure 6. In any cases, this density shouldn't be compared to the population density in the real world.

One should wonder why, after what was described in the former sections, we have chosen a so small space. There are mainly two answers, first because it was the limit of our best disposable screen capacity -for this first exploratory study it was necessary to keep an eye on the agent's behavior -, secondly because due to the overall huge memory requirement needed to simulate each simulation a thousand times, it was impossible to expect to be able to break this limit.

The density here play a fundamental role : as the agents change their displacement direction only when they encounter an other agent -which is obviously an oversimplified way to exchange information, which should be modified in a later version of the model -, if they are too few, an agent who has been assign a bad initial direction, that is a direction which doesn't put a firm on his trajectory, will have a very small probability to ever change direction and, by consequence, to find a job.

In revenge, as the population density rises, the number of direction change increases, offering a better chance to an agent to find a job. But every flow has its ebb. When the population density continue to rise, the agents can eventually change so often direction until to lengthen their search toward the failure to find a job. Theoretically, to every space dimension must an optimal population density be associated.

So, the way agents move and interact in their displacements is an important part of the model which condition the simulation results. This part of the model will be discussed later on this paper. No match can be fit outside an institutional arrangement, since institutional arrangements define the rules of the game. The first institutional arrangement preluding to the matches is the full-centralized one. The many second institutional arrangements are the partially-decentralized ones. The third institutional arrangements correspond to the markets.

The full centralized arrangement

In the full centralized arrangement there is one and only one unemployment agency which collects the labor demands and offers. This agency has a specific location randomly drawn in the interaction space which conditions the arrival of the labor demands, as well as the labor supplies, according to a specific distance -euclidian, Manhattan. . .calculated from the initial position of the firms and laborers.

According to the distances, two sets are sequentially constructed : the set of the demands and the set of the supplies. The constitution of those sets takes time, because of the arrival time of each demand as well as each supply but also because each agent's file must be studied by some public agent who, making seriously his job, uses a specific treatment to obtain the various A h and A i from the agents. For each set, if the information corresponding to two agents arrive at the same time they are ranked according to a random draw.Then two options are possible :

x the records are sequentially treated according to their place in the queue. If they matches, they get out it. If not, one of the two files is down graded of one rank and a new match is tried.

y in the second option, the public agent wait until a certain number of files of the two types are gathered. Then, he tries to find an optimal match among the two constituted files. If there is some agents that cannot be matched, they wait until the correct number of files is anew gathered. 21 In that case, there is no matter that the distance be euclidean or L 1 .

The partial centralized arrangement of degree n

One can observe that in many cases, there are not one labor agency but many, dispersed on the interaction space. This is because the central planer want to diminish the travel cost to the labor agency. Two arrangements are possible

x in the first arrangement, the agents are attracted by the labor agency which is the nearest of their location.

y in the second arrangement, each labor agency is dedicated to some tag intervals and agents must join a labor agency according to their A h and A i . One more time, one can calculate an index of efficiency an index of efficiency n O pc -more exactly two indexes, the first, n l O pc when the arrangement favors the location, the second, n t O pc when it favors the tags -, then one can construct the distances to the unaccessible O D -i.e.

: n l D f c = | n l O f c -O D | and n l D f c = | n t O f c -O D |.

The market arrangement

Of course, all this centrally and partially organized arrangements must be compared to the market one under highly imperfect information. In this simulation, the two sides of the market are differentially treated, in the sense that firms are fixed and only laborers move. Whichever be the interpretation of the interaction space, firms use to send messages to recruit not to move. In the contrary, the most part of the labor force uses to move to find a job. Of course, there are exception to those behavior where comportment are reversed but this is rare and more or less reserved to high level employment. That is not to say that firms adopt a passive mode of acting. They send messages informing of their characteristics the laborers who, in the curse of their moves, enter in their reception field which is raising simulation step after simulation step. This gave to firm an important advantage on laborers : to send a message corresponds to travel in all direction providing that messages diffuse uniformly in the interaction space. But here, the radius of the circle inside which the information of a firm is diffused, is finite and common to all firms. When receiving the information, a searching laborer can decide to change direction to reach the firm which has send the message if it announces a potentially better match than the one he was moving for, if he has yet receive a message -see figure 9.

Obviously, this can be a bad decision since the job can be already given to another laborer when he will reached the firm but it reflect the hazard of the search for a job.

Under the initialization stage, and whichever be the context of the simulation, the directions chosen by the laborers are drown at random in a uniform distribution. This direction is modified, again by a random draw, after any interaction with an other agent.

One must notice that in this paper, laborers can move freely in any direction all at the same displacement speed without to be restrained by any barriers or to be slowed in certain regions of the interaction space, restrictions which could be used to model some sociological weight that must disadvantage a certain part of the laborer population 22 . That is not to say that the modelling environment is not enough flexible to take in account those subtleties but here one concentrate the focus on the institutional organization of the matches. It's also the reason that their is no breaking of the realized 22 But we must keep this for a later paper.

matches. One want to abound in the case of a highly stylized market to the debate between centralization and free market taking in account the highly asymmetric level of lack of information from the agents.

An other point, on which it is necessary to insist, come from the space's size. Here, one poses it as fixed as it seems to be in the actual world. In reality, this is not really the case since the space can change over time through public work or many other operations.

Here, the space is more metaphoric than real. This has a direct consequence : when the agents' density increases, the number of collision grows accordingly. In the limit, for a very dense population of H-agents some could be prevented to move.

For each simulation, one must also construct an efficiency index -O m -, which in turn will be compared with the demiurge index -D D .

Results

To compute the results, one has to use a linear programming solver since one has shown earlier that the two sides matching problems are easily amenable to a linear program. The only problem was to make the link between the market simulation platform in such a way that each configuration can be saved an leads to the desired indicators.

First of all, as the characteristics of the population are randomly drawn -the two tags -, it was necessary to understand the simulation as some Monte-Carlo experiments: for each experiment, the scenario has been played 10000 times which gives a strong credibility to the results. One has successively worked with a couple ( H, I) equal to 10 × 10, 100 × 100,. . . , 1000 × 1000, 1200 × 1200 and finally 2000 × 2000. Of course, those dimensions are very small in comparison with the actual dimension of the true labor markets but, at least for those cases, all information is computable.

In all the cases, the protocol has began by the computation of the Pareto optimal matches, then one has used the greedy algorithm and one has run the moving agents platform.

In the case of all type of pairing, two results are of interest :

x the number of obtained pairs ;

y the level of the index of efficiency.

To evaluate the level of index efficiency, we will use the Hellinger statistic or more exactly the Hellinger distance between two discrete distributions p and q is define as :

H(p, q) = 1 √ 2 i ( √ p i - √ q i ) 2
which is no more, to the (2) -1/2 normalization constant23 , than the euclidian norm between the square root of the probability of realizations according to both distributions. Because of that it verifies the triangular inequality. One can show that :

H(p, q) = 1 - i √ p i q i
When for all i, p i = q i , H(p, q) = 0. It is equal to 1, whenever p i = 0 and q i = 0 and vice versa. So, smaller Hellinger distances give a better proximity between p and q. In what follows 4.0.1 The 10 × 10 market First of all, before anything else, it is not useless to call back that the optimum of Pareto doesn't look for the maximum number of matchings but the maximal value of the sum of surpluses. 10, that is the table obtained for 10 × 10 markets simulated 1000 times, it is obvious that for a little number of potential matches, the optimal matching linear programming algorithm -LP -largely dominate the two other ways to obtain an arrangement. Not only, in average it finds a greater number of matches, but it can find some perfect arrangements -all 10 jobs are finally allocated to all 10 workers -, when the greedy algorithm -GA -as well as the free search -MA -never find any one. That is to say, that when the two other ways to allocate the jobs are used, it is possible to stop on a non zero level of unemployment which should not be given the size of the market.

From the point of view of the total surpluses, not only it gives a greater mean than the two others, but one can also observed that when H-agents search by themselves the distribution of the surplus is left biased when it should be, according to the Paretooptimal matchings, right biased.

So, two important points should be noticed that :

x the greedy algorithm and the search are both highly inefficient even if the greedy algorithm -which is less demanding from the computational point of viewseems to be less inefficient.

y One can observe a very unrealistic level of unemployment -more than 40%. It seems to originate in the excessive requirement of the agents of which the mode of determination of the reservation tags creates an over-tension on the pairing possibilities. More specifically, this result seems to ensue from the fact that a lot of H-agents will agree to mate only with I-agents characterized by a greater tag that their own and vice versa.

One has two strategies to overcome this unwished effect:

x to change the value of ε h , μ h , ε i and μ i , mainly to distinguish them from each other and certainly to consider less demanding populations which always agree to mate with agents of tag lower than theirs -i.e. : ε h = 0, ε i = 0 and μ h = μ i = 0. y to change the distribution -until now, uniform -from which those parameters are drawn. One can for instance choose an exponential distribution for the H-agents and a gaussian one for the H-agents.

But one must understand that, that in both cases, the choices are highly arbitrary. That is not to say that it is uninteresting to try to do it. Now one can look at what happen when the dimension of the market rises.

4.0.2 From the 100 × 100 to the 2000 × 2000 market

In the annexe, one has collected some histograms showing the evolution of the matching problem according to the population size. What could be said from the histograms collected in the annexe is that from the point of view of the surpluses, even for small dimension the greedy algorithm seems to perform nicely. The most important part of the results can be gathered through the Hellinger distance in the figure 11. As the optimal matching is calculated on the surplus basis, the main results are linked to the value-Hellinger distance -upper graph in the figure 11. For what concerns the greedy algorithm, obviously, it is never able to implement on a 10 classes histogram -the value for which figure 11 has been implemented -the matching algorithm. But, it performs not so badly since in average the associated Hellinger distance is of 0.11 for a standard deviation of 0.073. That is to say that the greedy algorithm is biased, which is not an unexpected news, but the bias is more or less of constant dimension on the population size, which gives a measure of the performance between the first best and a second best optimum one could expect with it. This results is of interest: according to the complexity and the memory size needed to compute the results, it is clear that reaching the optimum is unattainable with any computer program if the population size is big enough. So, less complex and memory requiring, the greedy algorithm, even if one knows that it will also encounter some limits can be implemented to approach the optimum when it is not possible to compute it. More than that, for a given population size, the Hellinger distance often fluctuates very lightly with the number of classes used to compute the histograms and is nearly invariant for small size. But, when one exceed a size of 500 × 500, it nearly uniformly fall under 0.1. This shows a nearly adequation on the distribution of the optimal matches and the distribution of the greedy matches24 .

One can also observe that the simulation, that is letting agents searching and finding by themselves their jobs is a very bad solution, since its divergence from the optimum distribution always rises until reaching a complete divergence from it when the population size is 500. According to the fact that, in this model, agents don't really share any information local nether local nor global, it is not so surprising.

Looking at the second aspect of the computations, one can observe that the Hellinger distance for the simulated matches, with is initially falling to a minimum, finish to rise to the maximal divergence. This could be attributed to the invariance size of the space which finish, size by size, by being overpopulated. This leads to two much direction changes which, from a chance to increase the opportunities to meet companies, eventually represent an obstacle to the job search.

5 From one to many labor markets As underlined earlier, sooner or later one arrives to the limit of the optimum computation. In that case one solution is to break the problem in many subproblems. Formally this should create some problems since in dividing the space, one create artificial boundaries which are restructured as tores but which are not really. In that case, one has been astonished to see that the comparison between the 2000 × 2000 and 4 500 × 500 was extremely good with a Hellinger distance falling in a uniform way from 0.089442719 for a 10 classes histogram to 0 for a 40 classes histogram.

Conclusion

The main result of this study is that a population with a so much degree of requirement but managing so badly the information which it could have, should never be left to itself. In that case, the free market could not implement the Pareto optimum, contrary to the fact that lets believe the walrasian dominant approach.

A central planer should intervene and, by procuration, should manage the matching problem. Unfortunately, we have shown from the theoretical view point as well as the empirical view point that the program necessary to accomplish this task is, in all likelihood, non computable because of the size of the problem. Nevertheless, one can conjecture that, in breaking the main problem in subproblem, one can reach a second best optimum which is not so far from the true one. This is the solution to which one resorts in the real world. In this world, it seems not to be effective. In fact, the main difference between the real world and the simulated one is that contrary to the former, the later offer potentially one job to each demand. Obviously the matching of agents to tasks or job is not a simple task whenever their number is huge and one can be. To finish, one are conscious that the agents must be gifted with more capability in the treatment of the information 25 . This will be done in an other paper along with the many sociological barriers to reach employment.

Figure 1 :

 1 Figure 1: Tensions indicator on the French Labor Market -Quarterly Data

Figure 6 :

 6 Figure 6: Population density for a 500 × 800 grid

Finally, the center

  can calculate an index of efficiency O f c -in fact two indexes can be constructed O s f c for the case were the records are sequentially treated, and O op f c for the case where a partial optimum is attempted -which can be compared with the unaccessible O D in constructing a distance 21 D f c = |O f c -O D |.

Figure 7 :

 7 Figure 7: The Full Centralized Arrangement

Figure 8 :

 8 Figure 8: The partially Centralized Arrangement -n = 2

Figure 9 :

 9 Figure 9: Change of direction

Figure 11 :

 11 Figure 11: Hellinger distances to Optimum -for a 10 classes histogram

  It is exactly what, in France, Pôle Emploi does when it open many branches. For instance, one can look at what happen when 4 branches of 500 × 500 instead of one of 2000 × 2000, are opened.

Figure 12 :

 12 Figure 12: A Division of a Torus

An approximate Pareto optimality.

As an example, one can look at the literature on the Traveling Salesman Problem which shows that the optimal travel for a limited number of cities can put down the best computers programs -see[START_REF] Applegate | The Traveling Salesman Problem: A Computational Study[END_REF].

In 1940, Atanasoff -see[START_REF] Atanasoff | Computing machine for the solution of large systems of linear algebraic equations[END_REF] argued that it takes 8 hours to an expert to solve a system of 8 linear equations in 8 unknowns. It seems that the greatest system of equations ever solved by hand emerged on the occasion of the Principal Triangulation of Great Britain.[START_REF] Clark | Elements of Algebra[END_REF] notes a sub-system of 77 equations. Nowadays, computer scientists think that one cannot expect to solve a system of more than 10

if the matrix is full and perhaps 10 8 is there are exploitable special structures of the problem, which is not the case here.5 For instance the free Lp_Solve and GLPK or the commercial CPLEX, GUROBI or MOSEK which are freely distributed for academics.

It is not sure that this number includes the services.

The problem comes from the Uzawa equivalence theorem which states the equivalence between the general equilibrium and the Brouwer theorem -see[START_REF] Uzawa | Walras existence theorem and Brouwer's fixed point theorem[END_REF]. In fact,[START_REF] Hirsch | Exponential lower bounds for finding Brouwer fixed points[END_REF] have shown that the lower bound for worst-case computation of Brouwer fixed points is exponential in the dimension of the problem, that is the size of the commodity space in the Arrow-Debreu version of general economic equilibrium.

[START_REF] Anderson | Dynamic matching and evolving reputations[END_REF] demonstrate that the PAMs are not validated when dynamics are introduced. They specifically test PAMs in the presence of evolving public Bayesian reputations and stochastic production.

In this paper, one has retain an uniform drawing such as, for t h as well as t i , t {i}∨{h} U[8, 11].For the acceptation intervals one has set ε h = μ h = ε i = μ i = 0.1.

The demiurge is the deity responsible for the creation of the physical universe in diverse cosmogonies.

One has not been able to compute the optimal matching allocation with Lp_Solve for H+I = 1000 because the program has exhausted the memory in writing the transitory results.

This is a well known phenomenon revealed in the qaly -Quality Adjusted Life Years -assessment in health economics -see, among other[START_REF] Mehrez | Quality-adjusted life years, utility-theory, and healthy-years equivalents[END_REF].

Obviously,[START_REF] Sakoda | Minidoka: An Analysis of Changing Patterns of Social Interaction[END_REF] has been unnoticed before[START_REF] Schelling | Models of segregation[END_REF] and one can conjecture that it was a sentiment of frustration which decided the former to published on the subject, 22 years after his doctoral contribution.

The neighbourhood may be of the von Neumann style -i.e. : in crux to the first order neighbour cells -or more completely through the height direct neighbours.

To be more precise, in the visual version of the simulations, the continuity of the space is limited by the dimension of the pixel. But we have also a batch version for huge simulations in which the space is truly continuous, since there is dimension imposed support.

A priori, since one can imposes any type of restriction on the space.

In the representation of this space, one can keep a grid to give an artificial structure to the space which help to localize every agent.

It ensures that 0 ≤ H(p, q) ≤ 1.

One should notice that for all the other calculated distributions the Hellinger distances are stable and of high level.

Appendix

In this appendix, one first give the data needed to construct the figure 11, then one presents the histograms for the 100×100 market, the 500×500 market and the 1000×1000 market.

Data


Xulai Sun

Walter Kob

Raphael Blumenfeld

Hua Tong

Yujie Wang

Jie Zhang

Friction-controlled entropy-stability competition in granular systems

Using cyclic shear to drive a two dimensional granular system, we determine the structural characteristics for different inter-particle friction coefficients. These characteristics are the result of a competition between mechanical stability and entropy, with the latter's effect increasing with friction. We show that a parameter-free maximum-entropy argument alone predicts an exponential cell order distribution, with excellent agreement with the experimental observation. We show that friction only tunes the mean cell order and, consequently, the exponential decay rate and the packing fraction. We further show that cells, which can be very large in such systems, are short-lived, implying that our systems are liquid-like rather than glassy.

Dense granular materials show a highly complex response when subjected to repeated cycles of shear, mainly due to the strongly dissipative, hysteretic, and nonlinear interactions at the frictional contacts. The structures of such systems self-organize dynamically and show characteristics that on large scales appear to be universal [1]. Previous works focused on the motion on the particle scale [2,3] and phase behavior [4][5][6], the slow relaxation of stress and density [7][8][9], and the complex spatio-temporal dynamics [10][11][12]. These dynamics are relevant in a wide range of fields, including the aging and memory of glasses [13], fatigue of materials [START_REF] Suresh | Fatigue of materials[END_REF], catastrophic collapse of soils and sand in civil and geotechnical engineering [START_REF] Terzaghi | Soil mechanics in engineering practice[END_REF], as well as in geological processes, such as earthquakes and landslides [START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]. Despite this multitude of investigations, the nature and role of the contact network, a key quantity of granular systems, are far from being understood and in particular the relation between the properties of this network and the friction between the particles is currently not known. Intuitively, one expects higher friction to give rise to looser structures and, therefore, to overall larger cells in the contact network. It has been proposed that the characteristics of granular packs can be understood from entropic considerations [START_REF] Edwards | [END_REF]. In turn, the entropy should depend on the driving process and friction. To quantify this relation, we have created a set-up that allows to analyse entropy of contact networks for widely different friction coefficients.

We consider a two-dimensional (2D) system and focus on its cells, defined as the smallest (aka irreducible) closed loops of contacts between particles that touch at least two neighbors. Other than their importance for the mechanical properties of the system, the properties of cells affect heat conduction in granular assemblies and, in three dimensions, the permeability to fluid flow. Cells are closely related to quadrons [1,[18][19][20], which are the smallest volume elements of granular assemblies and play a fundamental role in granular statistical mechan- ics [START_REF] Edwards | [END_REF]18,21,22]. During cyclic shear, contacts are continually broken and made, leading to cells being created and annihilated. Here, we investigate the effect of friction on the structure and quasi-static dynamics. Our findings are the following: 1. The cell order distributions are exponential and can be derived from a maximum entropy argument. 2. The structural characteristics result from a competition between entropy and mechanical stability, which the friction tunes. Specifically, mechanical stability curbs entropy by increasingly limiting the number of possible cell configurations as the cell order increases, making entropy more dominant as the inter-particle friction increases. 3. Our systems are liquid-like rather than glassy.

The experimental setup, known as the stadium shear device [START_REF] Miller | AIP Conf Proc[END_REF], is sketched in Fig. 1(a). We used a stepping motor to drive periodically two stainless steel sprockets, connected to each other by a rubber belt which was corrugated on the inside to ensure no-slip between it and the particles. The particle assembly was sheared between the two parallel sections of the belt in the central region of the device and recirculate under the two sprockets. We observed the particles within the blue shaded area, shown in Fig. 1(a), which thus mimics simple shear between two infinite parallel boundaries [START_REF] Miller | AIP Conf Proc[END_REF]. We applied a cyclic strain, whose maximum varied from 3% to 10%. In comparison, the yielding strain of this system is around 3% [START_REF] Zheng | [END_REF]. To study the effects of friction, we used: gear-shaped nylon particles (friction coefficient µ → ∞) [25], a combination of stainless-steel cylindrical rings and gear particles, photoelastic disks (µ ≈ 0.7) [26], ABS plastic cylindrical rings (µ ≈ 0.32), stainless-steel cylindrical rings (µ ≈ 0.3), and Teflon-coated photoelastic disks (µ ≈ 0.15). To minimize crystallization, we used a 50%-50% binary mixture of particles of size ratio 1:1.4. The diameter of the small particles was 1.0 cm, except for gear particles, whose small particles had a pitch diameter of 1.6 cm and a tooth height of 0.36 cm. The total number of particles was about 2000 in each system. Each experimental run was started by depositing particles randomly inside the stadium and applying a cyclic shear until the system reached a steady state. Below we show that fewer than 100 cycles suffice to de-correlate the system from its initial state. Halting at every cycle the quasi-static process at maximum negative strain, we took a snapshot of the system, which allowed to monitor the static structure and the dynamics stroboscopicallly. At steady state, the resulting packing fractions were 0.74 ± 0.01 for the gears, 0.81 ± 0.01 for the photoelastic particles, and 0.83±0.01 for the Teflon-coated photoelastic, stainless steel, and plastic particles. These values are consistent with the expectation that in a sheared system at constant pressure, increasing the friction results in a decreasing packing fraction. More details on the experiment are given in the Supplementary Materials (SM) [27]. We show the results for a strain amplitude of 5% and the results for other values are qualitatively the same. For each given configuration of particles, we first identified the contacts and then constructed the cells, see Fig. 1(bc) [1,20]. Our results should be insensitive to the exact definition of cells [30,31]. Due to the finite resolution, we had to use a threshold to define inter-particle contacts. We have checked carefully that the results and conclusions presented here are not sensitive to this threshold, see SM [27] for details.

The cell order n is defined as the number of contacts that define the cell, see Fig. Figure 2 shows the cell order distribution (COD), p(n), for six systems with different friction coefficients. The graph demonstrates that p(n) is described very well by Cell Order Distributions of six different systems: Gears (triangles), gears/stainless steel cylindrical rings mixture (circles), photoelastic disks (squares), plastic cylindrical rings (inverted triangles), and Teflon-coated photoelastic disks (diamonds). The solid lines are the exponential functions predicted by Eq. ( 2). Inset: The variance of the COD for the different systems, with the same symbols as in the main panel and stars for photoelastic disks under pure shear. All measurements fall exactly on the analytical function predicted from the PDF in Eq. ( 2) The small spread of symbols within a given cluster of same-type particle represents variations of strain amplitudes from 3% to 10% and boundary pressure of ±30%.

an exponential function, p(n) ∝ e -n/y . We find that the exponential form is independent of the friction coefficient, maximum strain, particles stiffness, and whether we applied a simple or pure cyclic shear [START_REF] Zheng | [END_REF]. As expected, the higher the friction the larger the occurrence frequency of high-order cells -up to n = 30 in the gear system -enabling us to determine the distribution function to high accuracy. Our CODs are somewhat different from those found in numerical studies of granular systems, in which the sample was gradually tilted [31] or isotropically compressed [1,20], indicating that CODs: (i) are sensitive to driving protocol and (ii) differ between the quenched and quasi-static steady-states.

We can understand this result in the context of the granular statistical mechanics approach [START_REF] Edwards | [END_REF]18,21,32]. It is plausible that the structure self-organises via a competition between entropy, i.e. increasing disorder with the largest possible number of cell configurations, and a constraint of mechanical equilibrium, which excludes unstable configurations. Expecting higher friction to enable more high-order cell configurations, thus increasing the importance of entropy, we neglect the stability constraint in our systems and use a maximum-entropy argument to derive the steady-state COD, p(n). We impose two constraints: that p(n) is normalized and that it has a well-defined mean, n , represented by two Lagrangian multipliers, x and -1/y. Maximising the Gibbs entropy,

S = - ∞ n=3 p(n) log p(n) + x ∞ n=3 p(n) -1 - 1 y ∞ n=3 np(n) -n , (1)
yields

p(n) = (1 -e -1/y ) e -3/y e -n/y , (2)
with y = log n -2

n -3 -1
the typical decay of p(n). Using Euler's relation for planar graphs [START_REF] Coxeter | Regular Polytopes[END_REF] (see SM [27]), we also have

y = log 4 6-z -1
, with z > 2 the mean coordination number. To test this result, we determined for each system the value of n and used it to get the corresponding value of y. In Fig. 2, we include these theoretical predictions for p(n) from Eq. ( 2) and find that the PDFs fall on the predicted curves almost perfectlywithout any fitting parameter!

The independence of the theoretical PDF of higher n moments is non-trivial and reminiscent of the Boltzmann distribution, which depends only on the mean energy. In thermal statistical mechanics this is because the energy is defined up to an arbitrary constant [START_REF] Feynman | The Feynman lectures on physics[END_REF]. For the energy to be a proper extensive macro-quantity this constant must cancel on calculating its higher moments, which is only possible if the distribution is independent of the higher moments. In contrast, there is no physical reason why adding an arbitrary constant, n → n + n 0 , should not change the higher moments of p(n). To test whether or not higher moments need to be included via additional Lagrange multipliers, we compare the experimentally computed variance of the COD, for different strain amplitudes, to that calculated from Eq. (2). The inset of Fig. 2 shows that the two coincide perfectly, establishing that the variance of the COD depends only on n and higher moments need not be included.

In this derivation we also neglect spatial cell-order correlations and to test this assumption we measured the conditional probabilities that a cell of order n is neighbor of a cell of order m, p(n|m). In Fig. 3, we show these conditional PDFs for different values of m for gear particles, photoelastic disks, plastic rings, stainless-steel rings, and coated photoelastic disks. For all systems, the PDFs collapse onto a master curve, establishing that p (n|m) is independent of m and that hence the cell orders of neighboring cells are not correlated. This observation and the success of the maximum-entropy derivation support the view that local entropic effects play a key role in determining the cell structures.

A way to characterize cells geometry is by the distribution of their internal angles, θ. Figure 4 shows the conditional PDF of θ, given the cell order, q(θ|n), for three representative friction systems: the gear (the highest friction, µ → ∞), photoelastic (intermediate friction, µ = 0.7), and Teflon-coated photoelastic systems (low friction, µ = 0.15). The PDFs are similar for n ≤ 5, mainly due to geometric constraints: when n is small then, given the particle radii and their sequence around a cell, only n -3 internal angles are variable, limiting the explorable configuration space. In contrast, the PDF tails for θ > 180 o are increasingly suppressed with reducing friction. This corresponds to inhibiting too elongated cells [20], which suppresses disorder and further supports our thesis that friction tunes the entropy-stability competition. This competition between the roughly exponential increase of the total number of cell configurations with n and the restraining of too elongated cells also affects the distribution of cell aspect ratios. We studied this effect by investigating the ratio of the mean volume of n-cells to that of the equivalent regular n-polygon, γ(n) ≡ V c (n)/V rp (n). The more severe the inhibition of elongated cells the closer is γ to 1. In other words, increasing entropy reduces γ while the stability restrains this trend. Figure 5 demonstrates this competition and shows how friction modulates it: entropy dominates in high-friction systems, with γ falling monotonically, while in low-friction systems the stability constraint eventually forces γ up as n increase. The value of n, when γ starts increasing, depends on both friction and the driving protocol. For example, in the compression processes studied in [20] (also shown in Fig. 5), γ eventually increases for both low and high friction.

An important aspect of the dynamic structural selforganization is the lifetime of cells. Cells appear and disappear as contacts are made and broken [START_REF] Wanjura | [END_REF]. We de- fine S(t) -the probability that a cell existing at time t 0 neither merges with another cell nor splits until time t 0 + t, with time measured in units of cycles. As shown in Fig. 6, S(t) decays sharply, a trend consistent with the assumption that merging and splitting of cells is a local, uncorrelated process. Using exponential fits for guidance, the relaxation times are short: τ ≤ 0.84, ≤ 1.43 and ≤ 0.45 for the gear, photoelastic and Teflon-coated photoelastic particles, respectively. The cell lifetimes in the Teflon-coated photoelastic system are significantly shorter than in the high-friction systems, consistent with the increased mechanical fragility as inter-particle friction is lowered. These results support our observation that the steady state is reached within at most a few dozen of cycles. The decreasing lifetime with n is due to a decreasing mechanical stability and is another fingerprint of the interplay between process-governed and friction-tuned entropy-stability competition. Although the decay of S(t) is faster for small n in the photoelastic particles than in the gears, this trend is reversed for large n, consistent with the expected increased stability of high-order cells with friction. Including in Fig. 6 the survival probabilities of rattlers in these systems, we observe that they remain disconnected from the network for longer than the lifetime of any cell order, suggesting that, The ratio γ = Vc(n)/Vrp(n) increases to 1 as cells are more round and vice versa. Increasing entropy reduces γ, opposing the effect of mechanical stability. The monotonic decrease of γ(n) at high friction indicates that entropy dominates over stability, while the up-turn of γ(n) at low friction is evidence to the increasing restraining effect of stability. The latter behavior was also observed in [20] (included with permission), albeit under different dynamics.

to lowest order, their effect on the structural organization can be neglected in our dynamics. The rattlers survival time, while always longer than those of cells, also drops with decreasing friction indicating the rapid reconstruction and change of topology in low-friction systems. Our observations have another implication. A number of works in the literature model granular systems as glasses, especially in the quasi-static regime, where dy-namic processes are slow. Yet, the short relaxation times we observe suggest that, in our process, the medium behaves rather as a liquid. Whether this conclusion extends more generally to quasi-static dynamics at lower shear rates and how it depends on the density, especially near the jamming point, remain open questions. Intriguingly, our high-friction systems contain a non-negligible fraction of particle strands with z = 2 (see SM [27]). Similar strands, of dynamically cross-linked particles, have been observed in what is known as 'empty liquids' [36][37][38][39]. Whether this similarity provides a route to study empty liquids via macroscopic high-friction granular systems remains to be explored.

Our work shows that the value of the inter-particle friction coefficient does not affect the fundamental structural characteristics of quasi-static dynamic granular systems qualitatively. Rather, it modulates the entropystability competition, expanding the cell configurations space as it increases. This tunes the COD and affects the packing fraction and the cell survival time distributions. Finally, being able to understand the cell structural characteristics from entropy considerations supports the premise of Edwards and collaborators that granular systems can be described by entropy-based statistical mechanics [START_REF] Edwards | [END_REF]18,21,22]. XLS and JZ thank I. Procaccia for valuable discussions. XLS and JZ acknowledge the student innovation center at Shanghai Jiao Tong University. WK is member of the Insitut Universitaire de France. This work is supported by the NSFC (No.11774221 and No.11974238).

Experimental details

During the experiment, particles were cyclically sheared within the stadium shear device driven by a stepping motor. In our system, we defined strain as the ratio between the displacement of lateral boundary and the width of the shear device, i.e. the distance between the two lateral boundaries. The strain was varied from 3% to 10%. This maximum strain ensured that not too many particles would leave the camera window during a cycle, which could affect the statistics. The shearing was halted every cycle at maximum negative strain for two seconds and two snapshots were taken to improve the accuracy of the image analysis. The granular assembly was sheared by the shearing belt, which was kept under a constant confining pressure of 12.5N/m. We estimate the inertial number of the particle motion at O 10 -4 , allowing us to regard the process as quasi-static [28]. A snapshot of the system was taken every second to allow tracking of all the particle trajectories during each cycles.

Following the initial random placing of particles, several hundred shear cycles were applied; Fig. 6 of the main text shows that this is sufficiently long to decorrelate the system and to reach a steady state. Each data set presented in the main text and here was obtained in steady state from snapshots over approximately 4000 cycles in the gear and mixed particles systems, 400 cycles for the photoelastic and Teflon-coated photoelastic systems, 2000 cycles for the plastic particles system, and 1200 cycles for the stainless steel particles systems. The number of cycles in photoelastic and Teflon-coated photoelastic particles is smaller than the others because of a small, but finite, probability of the particles to buckle away from the plane, invalidating a large number of runs.

In the mixed particle systems, the number of each kind of particles (large or small, gear or steel) were identical. At steady state, the resulting packing fractions were 0.74 ± 0.01 for the gears, 0.81 ± 0.01 for the photoelastic particles, and 0.83±0.01 for the Teflon-coated photoelastic, stainless steel, and plastic particles. These values are consistent with the expectation that in a sheared system at constant pressure, increasing the friction results in a decreasing packing fraction.

Except for the gears, the inter-particle friction coefficients were measured as follows. Two identical arrays of same-size and same-material touching particles were aligned on flat bars and positioned vertically face to face, as shown in Fig. S1. An increasing horizontal force, measured by a commercial gauge, was applied to the top array until it started sliding when it reached F max . The friction coefficient was determined as the ratio of F max to the total weight of the top array.

Procedure to identify contacts

We used the center-to-center distance as the main parameter to detect a contact between two neighboring particles. To correct for possible lens distortion effects, we tiled every image into a grid of squares 80 pixels by 80 pixels. Within each tile of each image, we computed the local pair correlation function and then averaged it over images of the same tile in different cycles to obtain the ensemble-averaged local pair correlation function. The first and third peak of this function were then used to determine the local radii of the small and large particles, respectively.

If the center-to-center distance between neighboring particles fell below the sum of their local radii plus a tolerance value, the particles were deemed to be in contact. A gear local radius is its pitch radius. Once contacts are identified within one image, we checked that the contact positioning around each particle does not violate mechanical stability. For a given inter-particle friction coefficient, the friction angle at a contact between two particles is α f = tan -1 (µ), which is defined as the angle between the contact force and its normal component, as shown in Fig. S2. In order not to violate mechanical stability, the reflex angle specified between any two neighboring contact points on a given particle, as shown in Fig. S2, should not exceed α max = π + 2α f , otherwise at least one of the two neighboring contact point is a false contact, which leads to the direct identification of rattlers among particles of only two contact points. This process was repeated until no particles in the system violate mechanical stability. This check was particularly useful for particles with two or three contacts and allowed us to both eliminate all false contacts and remove potential rattlers. At our image resolution, the pitch radii of the large (small) gears were 43 (34) pixels and the radii of the large (small) photoelastic, Teflon-coated photoelastic, plastic, and stainless steel particles were 30 (22) pixels. For the results presented in this paper we set the tolerance values of 2 pixels for photoelastic and steel particles, 2.5 pixels for plastic particles, and 3 pixels for gear particles. To test the effect of the finite resolution on contact determination, we varied the tolerance between 2 and 4 pixels and checked the effect on the cell order distributions. We found that the exponential form, Fig. 2 in the main text, was unaffected and that the decay rate changed by at most 12%. As an example we show in Fig. S3 the case of the gear system, demonstrating that the value of the threshold does not influence the functional form of the distribution. We point out that the mentioned 12% change is an upper bound, at least for photoelastic and Teflon-coated photoelastic systems. For example, we used a series of snapshot from a pure shear experiment in [START_REF] Zheng | [END_REF] and applied the above contact detection algorithm. We found that the number of identified contact points is only 4% higher than the actual contacts identified using the combination of the above algorithm and the additional stress information from the force chain images.

Effect of particle-base friction

In principle, the friction forces between particles and the acrylic supporting base could affect the system dynamics and neglecting them requires justification. We measured the friction between each particle type and size and the base, µ b , and calculated the particle-base friction forces, f b = µ b mg, with m the particle's mass. Under a confining pressure P , the inter-particle forces are typically P D, with D the particle diameter. As can be seen from Table I, for a typical pressure, P = 12.5N/m, f b /P D is less than 0.06 for all particles and hence the particle-base I. Particle parameters and the ratio of particle-base friction to inter-particle forces. Note that P.E. represents the abbreviation for 'photoelastic'.

Since the particle-base friction depends on the particle masses, we had also to establish that it does not affect different-type particles sufficiently to take these into consideration. We therefore checked the mobility of the different-type particles in the system of gear/steel mixture. The mean squared displacement (MSD) for the different kind of particles are shown in Fig. S4. While the small steel particles are more mobile than the large ones at short times, the MSDs of the two types, and in fact of all four types, become identical after 30-40 cycles. Moreover, the MSD of the lighter gear particles is in fact smaller at short times than that of the steel particles, indicating that particle masses hardly affect their mobility. Thus, the long-time dynamics is independent of particle type and the particles-base friction can be neglected.

Effect of the stability-entropy competition

on cell volume

The structural properties of cells with a given order, such as shape and volume distributions, are governed by the competition between the mechanical stability and entropy. The shape of a cell can be quantified not only by the internal angles distribution, as done in the main text, Fig. 4, but also by the ratio of cell volumes to that of the equivalent regular n-polygon, Fig. 5 in the main text. In this section we thus discuss how the volume of a cell depends on its cell order n.

To start, we note that cell edges can be defined in two ways, both of which give the same cell identity and order, but different volumes.

Definition 1: The edges extend between neighboring particle centers around the smallest (aka irreducible) void loops, e.g. the brown regions in S5(a). This definition useful mainly for disks and particles with aspect ratio to one.

Definition 2: The edges extend between contact points around such loops, blue regions in S5(b). This definition yields smaller cells and it has two main advantages: it is not sensitive to the particle shape and it describes directly the contact network. Since the latter underlies stress transmission during the dynamics of granular materials, as well as in the static states of granular materials, definition 2 has been shown to be useful for the analysis of granular structures [18,20] in that it allows one to define the smallest possible volume element in granular systems -the 'quadrons' [18]. A quadron is based on the contact network and both a particle and one of its neighboring cells contribute to its volume. To construct a quadron (see also Fig. S5(b)): (i) determine the centroid of a particle as the mean position of its contacts, ρ g ; (ii) determine the centroid of one of its neighboring cells as the mean position of the contacts around Photoelastic (μ=0.7) Gear (μ→∞) it, ρ c ; (iii) draw a vector between the two contacts of the particle that are shared with this cell, r cg ; (iv) the quadron is the quadrilateral, whose diagonals are r cg and ρ c -ρ g . An example of a quadron, shaded grey, is shown in Fig. S5(b). The quadrons tessellate the system perfectly, except for wildly non-convex particles, in which case a small percentage of them may overlap [29].

To analyze the cell volumes, we define the following quantities: V c , the volume of a specific cell as given by definition 2; V q , the sum over the volumes of the quadrons 'belonging' to a specific cell; V T , the sum over the vol-umes of the triangles, each of which is defined by the centroid of a particle around the cell and the two contact points that the particles shares with it, or in other words, the part of the quadron within the particle (see Fig. S5(b). By these definitions, V q = V c + V T for each cell.

To test the dependence of the volumes on the cell order n, we plot the cell volume per contact, V c /n, and the mean quadron volumes, V q /n against n for the gear (Fig. S5(c)) and photoelastic (Fig. S5(d)) systems. In these plots, volumes are measured in units of π R2 , with R ≡ (R s + R l ) /2, with R s and R l the radii of the small and large particles, respectively. The mean of V T increases linearly with n, see insets, which is to be expected with this distribution of particle sizes.

For these high friction systems an apparently reasonable fit for V c (n) is

V c (n) = n[A ln(n) -B] , (3)
with A and B constants, although we know of no theoretical model that predicts this form. The trend in high-friction systems shown in Figs. S5 can be understood as follows: At low values of n, geometry and mechanical stability constrain the mean volume of cells to approximate the equivalent regular polygons, V c (n) ≈ V rp (n) = n R2 cos 3 (π/n)/ sin(π/n), with some stretching along shear-driven tilted principal stress.

As n increases, more elongated cells can be realized, which suppresses the volume rate.

The increase in the fraction elongated cells at given cell order n with inter-particle friction leads to a decrease in the average cell volume, an effect that is seen clearly in Fig. S4: the higher the friction the slower the increase of the average cell volume with n. For example, at n = 7 the value of v c /n is just under 0.6 for the photoelastic system (µ = 0.7) while it is below 0.4 for the gears (µ → ∞).

In Fig. S7, we plot V c (n) for all the systems considered. As expected, reducing the intergranular friction steepens the curves at large n and they tend toward the limit of the regular polygon, included in the graph as well. While V rp ∼ n 2 for large n, we observe that, intriguingly, V c increases almost linearly with n for high-friction systems. This is again evidence for the ubiquity of elongated cells in these systems, which are suppressed by the mechanical stability constraint in low-friction systems. This rapid increase of the configurations space in high-friction systems, in spite of the stability constraint, is the reason for the domination of entropy in determining the cell structures distribution.

Euler relation and boundary corrections

To check for potential influence of the boundary on our results, we first use Euler's relation to establish a relation between the particles mean coordination number, z , and the mean cell order, n

n = 2 z z -2 + O N -1/2 . ( 4
)
Here N is the total number of particles and the second term on the right hand side represents boundary corrections. We then plot in Fig. S8 relation (4) for N → ∞ (solid line) and our measured values of z and n (symbols). Given sufficient statistics, which we have, any deviation between the experimental data and the curve would be evidence for existing finite boundary corrections. The experimental data falls squarely on the theoretical curve, showing that the boundary corrections to our observations are negligibly small. 4). The symbols represent measurements of z and n , corresponding to different systems (see legend) in different conditions: strain amplitudes, etc. The solid line is Euler's relation for N → ∞, in which limit the boundary correction vanishes. We observe that the experimental data fall nicely on top of the theoretical limit, establishing that the boundary effects are negligible. The small spread of symbols within a given cluster of same-type particle represents variations of strain amplitudes from 3% to 10% and boundary pressure of ±30%. for different slices of the gear system at a strain amplitude of 5%. The symbol colors correspond to the slice colors in the inset. No significant dependence of the PDFs on their position is observed. The solid black line is the average PDF of the system.

Homogeneity of the system

In a setup for simple shear there is a danger of a density gradient developing between the boundaries. To check the structure for statistical homogeneity, we divided the system into nine slices parallel to the shearing boundaries, as shown in the inset of Fig. S9, and calculated the COD in each slice. Plotting all the PDFs in Fig. S9, we observe that these are all exponential for n > 5, p(n) ∼ e -n/y , with y depending very weakly on the slice's position. To a good accuracy, we can thus regard the PDFs as collapsing onto a master curve. The noise for very large values of n we attribute to the decreasing statistics in this range of n. This collapse is consistent with the system being statistically homogeneous.

Internal angle distributions

Similar to the case of Teflon coated photoelastic disks show in the main text, Fig. 4, the internal angle distributions in the other two low-friction systems, plastic rings (µ = 0.32) and stainless steel rings (µ = 0.3), are qualitatively similar (see Fig. S10). As is the case for the low-friction systems discussed in the main text, the PDF tail at angles larger than 180 • is strongly suppressed due to the reduced mechanical stability of such cells. This is consistent with the plot in Fig. 5 in the main text, which shows that large cells in low-friction systems are more rounded.

Cells structural properties

A key characteristic of the local structure is the distribution of the particles' coordination numbers, C(z). Since these distributions are expected to depend on the parti- cle sizes, we show in Fig. S11 the C(z) for the large and small particles separately. In the gear system the mean coordination numbers for the small and large particles are z s = 3.09±0.07 and z l = 3.51±0.09, respectively, Fig. S11(a). These values are significantly lower than the densest state, z = 6, and are close to the marginally rigid limit of an infinite system, which indicates a very dilute structure. Figure S11(a) also shows that for about 10% of the particles z ≤ 1these are rattlers that do not transmit forces. The corresponding mean coordination numbers for the photoelastic particles are z s = 4.17 ± 0.03 and z l = 4.79 ± 0.03, respectively. These systems are denser than the gear systems, which is to be expected because of the relatively lower inter-particle friction. Lowering the friction further leads to smaller cells and increases the density and the mean coordination numbers: For the plastic particles we find z s = 4.20 ± 0.02 and z l = 4.79 ± 0.02, for the steel particles we have z s = 4.21 ± 0.03 and z l = 4.78 ± 0.03, and for the Teflon-coated photoelastic particles we have z s = 4.17 ± 0.03 and z l = 4.87 ± 0.03.

The fast structural relaxation of the system discussed in the main text suggests that these systems are liquidlike rather than glassy and because of the low density they resemble an empty liquids [36][37][38][39]. A key feature of empty liquids is their internal structure, which consists of strands of coordination number z = 2, whose lengths have been found to be given by an exponential distribution. To check for this feature in our systems, we have calculated q(l), the probability that a strand has exactly length l, as sketched in Fig. S12(a). In Fig. S12(b) we show that, as in empty liquids, q(l) is an exponential: q(l) ∝ exp(-l/l 0 ). For all our systems we find that the decay length is small, l 0 = 0.402, 0.156, 0.153, 0.145 and 0.139 for our gear, photoelastic, Teflon-coated photoelastic, steel and plastic systems, respectively, thus much shorter than the values found in empty liquids strands [36][37][38][39]. However, since l 0 depends on the details of the experimental setup, such as particle-base and interpaticle frictions, confining pressure, and possibly others, we cannot rule out the possibility that longer strands and lower densities may be achievable in other setups, which could then be used as convenient macroscopic models for empty liquids. 9. Cell survival probabilities of plastic ring and stainless steel ring systems Here, we show that the survival probabilities of cell orders follow the same pattern as for the systems of gear, photoelastic, and coated photoelastic particle systems, discussed in the main text. In Fig. S13 we show the probability that a cell present at time t = 0 has neither split nor merged with other cells at time t. These probabilities in the plastic and stainless-steel ring systems are qualitatively similar to the case of the Teflon-coated photoelastic particles shown in Fig. 6 in the main text. Note that the lifetime of rattlers in the steel particle system is significantly larger than that in plastic particle systems, indicating that an increase in the friction between the particles and the supporting base may slow down the integration of rattlers into the force carrying network.

FIG. 1 .

 1 FIG. 1. (a)The stadium shear device is mounted on a horizontal glass plate. Particles occupy the region within the stadium, including areas under the sprockets. A constant boundary confining pressure is maintained by two aluminium bars, mounted on four air bearings with constant forces pushing against the lateral walls of the belt. The imaging system includes a CCD camera and a LED panel above and below the device, respectively. (b) A cell is defined as the smallest loop of contacts. Two cells are drawn in the diagram. (c) A part of the gear system, containing cells, shaded blue..

  1(b) and (c) and the SM [27].

  FIG. 2.Cell Order Distributions of six different systems: Gears (triangles), gears/stainless steel cylindrical rings mixture (circles), photoelastic disks (squares), plastic cylindrical rings (inverted triangles), and Teflon-coated photoelastic disks (diamonds). The solid lines are the exponential functions predicted by Eq. (2). Inset: The variance of the COD for the different systems, with the same symbols as in the main panel and stars for photoelastic disks under pure shear. All measurements fall exactly on the analytical function predicted from the PDF in Eq. (2) The small spread of symbols within a given cluster of same-type particle represents variations of strain amplitudes from 3% to 10% and boundary pressure of ±30%.

2 FIG. 3 .

 23 FIG. 3. The conditional CODs p (n|m) for different systems. For clarity, data of stainless-steel particles, plastic particles, photoelastic particles and gears have been shifted vertically by factors of 5, 25, 125, and 625, respectively.

1 FIG. 4 .

 14 Photoelastic(μ=0.7)

  FIG. 5.The ratio γ = Vc(n)/Vrp(n) increases to 1 as cells are more round and vice versa. Increasing entropy reduces γ, opposing the effect of mechanical stability. The monotonic decrease of γ(n) at high friction indicates that entropy dominates over stability, while the up-turn of γ(n) at low friction is evidence to the increasing restraining effect of stability. The latter behavior was also observed in[20] (included with permission), albeit under different dynamics.

FIG. 6 .

 6 FIG. 6. (a-c) The survival probability of cells of different orders, n, and rattlers, in the systems of (a) gears, (b) photoelastic particles, and (c) coated photoelastic particles. (d) The corresponding characteristic relaxation times, τ (n), estimated from an exponential fit of the data in the main panel.

Figure S1 .

 S1 Figure S1. Sketch of the set-up for measuring the interparticle friction coefficients.

Figure S2 .

 S2 Figure S2. Sketch of a two-contact particle in a critical forcebalance state.

Figure S3 .

 S3 Figure S3. Effect of the finite resolution of contact determination on the cell statistics in the gear system. Cell order distributions calculated using different error tolerance values of 2, 3, and 4 pixels as specified in legend, show no qualitative difference, with the decay rate changed by at most 12%.

Figure S4 .

 S4 Figure S4. Mean squared displacement of the four different types of particles in the mixed system. The dashed line is a power law with slope 1.0.

Figure

  FigureS5. (a) Cells (shaded brown) defined by connecting particle centers. (b) Cells (shaded blue) defined by connecting particle contacts. The grey shaded region is an example of a quadron -the smallest volume element of the structure[18]. (c) Cell and quadron volumes for the gear system: The main panel displays the average volumes of cells per cell edge (circles) and of quadrons (square) as a function of n with a logarithmic scale. Inset: n-dependence (on a linear scale) of the volumes in the main panel, multiplied by n, and the the difference between them, VT (n) (triangles). (d) As in (c) for the photoelastic particles.

Figure S6 .

 S6 Figure S6. Cell and quadron volumes for the low-friction systems. The main panel displays the average volumes of cells per cell edge (circles) and of quadrons (squares) as a function of n with a logarithmic scale. Inset: n-dependence (on a linear scale) of the volumes in the main panel, multiplied by n, and the difference between them, VT (n) (triangles). (a) Plastic particles. (b) Steel particles. (c) Teflon-coated photoelastic particles.

Figure S7 .

 S7 Figure S7. The n-dependence of the mean cell volume for the various experimental systems and of the volume of the regular n-polygon. The volumes are measured in units of π R2 .

Figure S8 .

 S8 Figure S8. Check of the boundary correction, using Euler's relation, Eq. (4). The symbols represent measurements of z and n , corresponding to different systems (see legend) in different conditions: strain amplitudes, etc. The solid line is Euler's relation for N → ∞, in which limit the boundary correction vanishes. We observe that the experimental data fall nicely on top of the theoretical limit, establishing that the boundary effects are negligible. The small spread of symbols within a given cluster of same-type particle represents variations of strain amplitudes from 3% to 10% and boundary pressure of ±30%.

Figure S9 .

 S9 Figure S9. Cell order probability density functions (PDFs) for different slices of the gear system at a strain amplitude of 5%. The symbol colors correspond to the slice colors in the inset. No significant dependence of the PDFs on their position is observed. The solid black line is the average PDF of the system.

  Photoelastic(μ=0.7)

Figure S10 .

 S10 Figure S10. Distributions of the internal angle for different cell orders n for various systems. (a) Gear system, (b) Photoelastic system, (c) Plastic particle system, (d) Steel particle system, (e) Teflon-coated photoelastic system.

  Figure S11. Distribution of the coordination numbers for the small and large particles. (a) gear, (b) photoelastic (PE) particles, (c) plastic ring systems, (d) steel particles, and (e) Teflon-coated photoelastic particles. The shear strain is 5%.

Figure

  Figure S12. (a) Sketch of a strand of length l. (b) The strand length PDFs for the different systems. The strain amplitude is 5%.

  Figure S13.The survival probability of cells of different orders n and rattlers in (a) Plastic particle system and (b) Steel particle system.

Supplementary Material:

In this Supplementary material we provide additional information on: details of the experimental setup; the procedure to identify contacts; effect of particle-support friction; the cell volume calculation; Euler relation; the effects of boundary conditions on the structure; the cell internal angle distributions and cell survival probabilities in the systems of plastic rings and stainless steel rings; and the structural properties of the system.