text:
# Disappearance of Transverse Flow in Central Collisions for Heavier Nuclei

## Abstract

Using the Quantum Molecular Dynamics model, we analyze the disappearance of flow in heavier colliding nuclei. A power law mass dependence ($`\frac{1}{\sqrt{A}}`$) is obtained in all cases. Our results are in excellent agreement with experimental data which allows us to predict the balance energy for 238U + 238U collision around 37-39 MeV/nucleon.

Thirty years ago it was predicted by Scheid and Greiner Sch68  that in heavy ion reactions the nuclei will be compressed and heated and that this yields for non central reactions to in-plane flow $`p_x^{dir}`$. More than a decade later, this conjecture was confirmed by the Plastic Ball group pb . In the following investigations it turned out that this in-plane flow carries information on the nuclear equation of state ho . If the nuclear equation of state is stiffer, more compressional energy will be stored in semi-central reactions and, when released, more in-plane flow will be given to the nucleons.

The maximal density which is reached in a reaction depends on the beam energy as well as on the system size. The lower the beam energy the less is the compression. At very low energies, the repulsive part of the nuclear equation of state, which appears at densities above the normal nuclear matter density, is not tested anymore and the nucleons feel only the attractive mean field. A typical example is the deep inelastic reactions in which the two nuclei rotate around a common center. This rotation creates in-plane flow as well but in opposite direction: Due to the common rotation the nucleons stick together for a while and will be emitted into the direction opposite to the impact parameter whereas the in-plane flow which is caused by compression will be in the direction of the impact parameter.

There is a beam energy at which the in-plane flow disappears when changing from the direction into that opposite to the impact parameter. It has been shown in the simulation of heavy ion reactions that this beam energy called balance energy, $`E_{bal}`$, 3a ; 1a ; 2a  depends on the nucleon-nucleon (nn) cross-section in the medium as well as on the potential 1a ; 2a . With the very recently measured $`E_{bal}`$ in 197Au + 197Au collisions 4a  (earlier only estimated values were available 5a ), a renewed interest has emerged in the field 6a . In addition to the Au system, balance energies $`E_{bal}`$ of 12C + 12C 7a , 20Ne + 27Al 7a , 36Ar + 27Al 8a , 40Ar + 27Al 9a , 40Ar + 45Sc 4a ; 7a ; 14a , 40Ar + 51V 2a ; 10a , 64Zn + 27Al 11a , 40Ar + 58Ni 6a , 64Zn + 48Ti 8a , 58Ni + 58Ni 4a ; 6a , 64Zn + 58Ni 8a , 86Kr + 93Nb 4a ; 7a , 93Nb + 93Nb 3a , 129Xe + 118Sn 6a  and 139La + 139La 3a  are also available. It is worth mentioning that most of the above studies were for the central collisions only. A few, however, also searched for the impact parameter dependence of the balance energy 4a ; 9a ; 14a ; 11a .

Apart from the directed in-plane flow, differential as well as elliptic flow has also been predicted very recently 12a .

The measurements of the balance energy over wide range of system sizes provide an excellent opportunity to pin down the role of the mass dependence, where only preliminary studies 4a ; 7a  have been performed yet. These preliminary studies suggest a power law dependence $`A^\tau `$ of the balance energy on the mass number of the system. Interestingly, most of the theoretical studies are done within the Boltzmann-Uehling-Uhlenbeck (BUU) model 3a ; 1a ; 2a ; 4a ; 7a ; 9a ; 11a ; 12a ; 13a ; 15a ; 16a ; 17a ; 18a . Some attempts, however, also exist within the framework of Quantum Molecular Dynamics (QMD) model 14a ; 19a ; soff95 ; 21a . Heavy systems are rather rarely analyzed in these approaches.

Our present aim is therefore to study the mass dependence of the balance energy in heavy colliding nuclei and to predict for the first time the disappearance of the collective in-plane flow in central 238U + 238U collision. We shall show that the mass dependence of $`E_{bal}`$ for heavier nuclei scales approximately more as $`\frac{1}{\sqrt{A}}`$ rather than as $`A^{\frac{1}{3}}`$ as has been suggested for light and medium colliding nuclei 7a . The present study is made within the framework of QMD model which is described in detail in refs. 19a ; soff95 ; 21a ; 22a ; 23a .

In the QMD model, each nucleon propagates under the influence of mutual interactions. The propagation is governed by the classical equations of motion:

$$\dot{𝐫}_i=\frac{H}{𝐩_i};\dot{𝐩}_i=\frac{H}{𝐫_i},$$
(1)
where H stands for the Hamiltonian which is given by:

$$H=\underset{i}{\overset{A}{}}\frac{𝐩_i^2}{2m_i}+\underset{i}{\overset{A}{}}(V_i^{Skyrme}+V_i^{Yuk}+V_i^{Coul}).$$
(2)
Here $`V_i^{Skyrme}`$, $`V_i^{Yuk}`$ and $`V_i^{Coul}`$ are, respectively, the Skyrme, Yukawa and Coulomb potentials. Momentum dependent interactions are not important at these low beam energies as Zhou et al. 18a  have shown. It is worth mentioning that both the soft and hard equations of state have been employed in the literature. Following 3a ; 1a ; 9a ; 14a ; 11a ; 12a ; 17a ; 18a ; 19a ; soff95 ; 21a ; 24a , we shall also use a hard equation of state. For the nn cross-section we chose an isotropic and energy independent cross-section of 40 mb for the present analysis. This seems to us a reasonable choice in view of the fact that cross-sections based on G-Matrix calculations differ widely and because at this energy most of the collisions are Pauli blocked. Hence different cross-sections do not produce a large effect if they are not too small.

Since we plan to study the heavier colliding nuclei, different nn cross-sections should not have much effect 4a . Further, it has been shown in ref. 15a  that the nucleons in the present energy domain ($``$ 80 MeV/nucleon) collide with average $`\sigma `$=55 mb. The isotropy of the cross-section also does not affect the reaction dynamics 25a . Similar assumptions were also made in refs. 1a ; 9a ; 11a ; 17a ; 21a ; 26a ; 27a . One should, however, keep in the mind that different nn cross-sections may affect the dynamics in lighter systems.

Using the above description, we simulated the central collisions with system mass A ($`=A_T+A_P`$; $`A_T`$ being the target mass, and $`A_P`$ being the projectile mass) $`175`$. In particular, we simulated 86Kr + 93Nb (b = 4.07 fm) 7a , 93Nb + 93Nb (b = 3.104 fm) 3a , 129Xe + 118Sn (b = 0-3 fm) 6a , 139La + 139La (b = 3.549 fm) 3a , 197Au + 197Au (b = 2.5 fm) 4a  and 238U + 238U (b = 0-3 fm) at incident energies between 30 MeV/nucleon and 80 MeV/nucleon at a step of 10 MeV. A straight line interpolation between steps was used to calculate the energy of vanishing flow $`E_{bal}`$. The reaction was followed till transverse flow saturates which is close to 300 fm/c for heavier colliding nuclei whereas it is $`200`$ fm/c for lighter colliding nuclei.

In fig. 1, we display the average directed transverse momentum $`p_x^{dir}`$ defined as:

$$p_x^{dir}=\frac{1}{A}\underset{i=1}{\overset{A}{}}sign\{Y(i)\}p_x(i),$$
(3)
where $`Y(i)`$ and $`p_x(i)`$ are, respectively, the rapidity and the transverse momentum of the $`\mathrm{𝑖𝑡ℎ}`$ particle for the reactions 93Nb + 93Nb, 139La + 139La, 197Au + 197Au and 238U + 238U. We see that the collective in-plane flow for 93Nb + 93Nb changes sign between 55-60 MeV/nucleon whereas it is already positive around 40 MeV/nucleon for 238U + 238U system. Further, the saturation time increases with the size of the colliding nuclei. The early onset of the flow in heavier colliding nuclei is due to the Coulomb forces that are much stronger in heavier systems compared to lighter nuclei. In addition, a large collision rate in heavier systems also contributes towards the early onset of the flow.

In fig. 2, we display the energy of vanishing flow ($`E_{bal}`$) as a function of the combined mass of the system. As stated earlier, $`E_{bal}`$ is extracted using a straight line interpolation between the calculated values of the in-plane flow. Here open squares represent our calculations whereas solid stars are the experimental findings. The dotted and dash-double-dotted lines are, respectively, the power law fits ($`A^\tau `$) to the theoretical values including 238U + 238U in one case and excluding 238U + 238U in the other case. The fit to the experimental points is represented by a solid line. All fits are obtained with $`\chi ^2`$ minimization. The values of $`E_{bal}`$ (obtained with a stiff equation of state and 40 mb cross-section) are very close to the experimentally measured $`E_{bal}`$. The fit to the experimental data yields $`\tau =0.52451\pm 0.06261`$, whereas that to theory yields $`\tau =0.53261\pm 0.16373`$. Once 238U + 238U is included, the $`\tau `$ decreases to $`0.50872\pm 0.10883`$. In other words, we observe a $`\frac{1}{\sqrt{A}}`$ dependence in the $`E_{bal}`$. Similar dependence can also be obtained with a least square fit. Based on the above findings, we predict $`E_{bal}`$ for the central 238U + 238U reaction around 37-39 MeV/nucleon (according to power law fits, it is around 37-38 MeV/nucleon whereas QMD simulation predicts around 39 MeV/nucleon). It is worth mentioning that most of the earlier mass dependence calculations 4a ; 7a ; 18a  could not reproduce the experimentally extracted slopes 4a ; 7a . Zhou et al. 18a  could reproduce the slope, however their analysis was only done for lighter nuclei $`200`$. Our calculations can reproduce the experimentally extracted slope very closely, therefore, we can predict the energy of vanishing flow in the 238U + 238U system. The large deviation of our $`\tau `$ value from the standard value ($`\frac{1}{3}`$) reflects the increasing importance of the Coulomb repulsion with the size of the system, as noted in ref. 4a . There the $`\tau `$ value is close to $`\frac{1}{3}`$ for masses $`200`$ 7a ; 18a , whereas it increases to $`0.45`$ when heavier systems like 139La + 139La and 197Au + 197Au are included. In the present analysis, we took only heavier nuclei ($`A175`$). Therefore the slope is steeper than the above cited values. If one also takes the lighter nuclei into consideration, our value also decreases to $`0.4`$. In other words, for lighter and medium nuclei, the balance energy $`E_{bal}`$ emerges due to the interplay between the mean field and nucleon scattering. However, for heavier colliding nuclei, Coulomb interaction is as well an important factor. It is of interest to see the contributions of the mean field (that includes the Coulomb interaction) and nn collisions towards the transverse flow at the balance energy $`E_{bal}`$. Following ref. 19a , we decomposed both these contributions in the simulations itself. At each time step during the reaction, the momentum transferred due to the two- body collisions and mutual mean field potential is calculated separately using equation (3). The separation was done at each simulated incident energy and a straight line interpolation was used. The decomposition is plotted in fig. 3 as a function of the total mass of the system. We find that the contribution of mean field towards transverse momentum is negative whereas it turns repulsive for the collision part. Both these contributions again obey a power law behavior with $`\tau =0.70931\pm 0.24234`$.

Summarizing, we present the disappearance of flow in heavier colliding nuclei with a prediction of balance energy for 238U + 238U around 37-39 MeV/nucleon. Our calculations (with a stiff equation of state and $`\sigma =40`$ mb) are in a very close agreement with experimentally extracted values ($`\tau _{th}=0.53261\pm 0.16373`$; $`\tau _{expt}=0.52451\pm 0.06261`$). Both these findings suggest a power law mass dependence $`\frac{1}{\sqrt{A}}`$. The contribution of the mean field towards flow at $`E_{bal}`$ is negative whereas it is positive for the collision part. Both contributions can be parameterized in terms of a power law.

This work is supported by the grant (No. SP/S2/K-21/96) from the Department of Science and Technology, Government of India.

Figure Captions

FIG. 1. The time evolution of $`p_x^{dir}`$ for four different reactions: (a) 93Nb + 93Nb, (b) 139La + 139La, (c) 197Au + 197Au and (d) 238U + 238U. Here a stiff equation of state along with constant nn cross-section of 40 mb strength is used.

FIG. 2. The $`E_{bal}`$ as a function of the total mass of the system. Solid stars are the experimental data whereas open squares are the present theoretical results. The solid line is a $`\chi ^2`$ minimization fit of power law ($`A^\tau `$) for experimental data whereas dash-double-dotted line is a fit for the corresponding theoretical result. The theoretical fit that includes 238U + 238U reaction is represented by dotted line.

FIG. 3. The decomposition of $`p_x^{dir}`$ at $`E_{bal}`$ into mutual mean field part and collision part as a function of the system size. The lines are the $`\chi ^2`$ fit of power law $`c.A^\tau `$.


format: text

text:
# Multi-partite Quantum Entanglement versus Randomization: Fair and Unbiased Leader Election in Networks

## 1 Introduction

Quantum entanglement is one of the most fascinating aspects of quantum physics. Two particles in an entangled state behave as if possesing a common state, even if they are physically separated over a great distance. Two maximally entangled qubits can be in the Bell state $`|\mathrm{\Phi }^+>=\frac{1}{\sqrt{2}}(|00>+|11>)`$. Each of the qubits here is in a mixed state with probability $`\frac{1}{2}`$ in state $`|0>`$ and with probability $`\frac{1}{2}`$ in state $`|1>`$. However, both together form a pure state . Such a pair of qubits is called an EPR pair (after Einstein, Podolsky and Rosen) . Such entangled qubits can exhibit what physicists call nonlocal effects. Nonlocal effects cannot occur between particles as per laws of classical physics unless communication is permitted between the physically separated particles, in which case communication would have to be at a speed exceeding that of light. Such behaviour was referred to as spooky actions at a distance in  and was earlier mentioned in the celebrated 1935 paper by Einstein, Podolsky and Rosen . Nonlocal effects however do not involve any communication and therefore the question of faster than light communication does not arise. Even if two agents Alice and Bob share entangled pairs, there is no way for Alice to manipulate her own particles to convey any information, she has, to Bob; information must be sent explicitly from Alice to Bob. In this context, it is necessary to investigate whether quantum computing resources help Alice and Bob solve a problem together with lesser communication complexity than that required when no quantum resource is used.

In Yao’s model , Alice and Bob are allowed to communicate using qubits. Kremer’s result in  does not show superiority of quantum communication over classical communication for the INNER PRODUCT function $`f(x,y)=x.y`$ where $`x`$ and $`y`$ are $`n`$-bit vectors each, $`x`$ given only to Alice as input and $`y`$ given only to Bob as input; the lower bound on qubits required for computing INNER PRODUCT is $`\mathrm{\Omega }(n)`$. However, the result of Buhrman, Cleve and Wigderson  requires only $`O(\sqrt{n}\mathrm{log}n)`$ qubit communication cost for the DISJOINT function, whose classical probabilistic communication complexity is $`\mathrm{\Omega }(n)`$ . The DISJOINT(x,y) function gives 0 if any bit of $`n`$-bit vector $`x`$ is equal to the corresponding bit of the $`n`$-bit vector $`y`$; otherwise, it returns a 1.

Just as two particles at a distance could be quantum entangled forming an EPR pair, it is also possible to entangle multiple mutually physically separated particles. One example is due to Greenberger, Horne and Zeilinger ; here, three particles are quantum entangled. Singh, Kumar and Pal  develop two protocols that entangle three nodes $`A,B`$ and $`C`$ in a network in the GHZ state $`(|0_A0_B0_C>+|1_A1_B1_C>)/\sqrt{2}`$ if two pairs (say $`(A,B)`$ and $`(B,C)`$) were earlier entangled as two EPR pairs. The first protocol uses the standard teleportation circuit  and the second one uses a new circuit developed in . Both protocols use only classical communication. Buhrman, Cleve and van Dam , make use of three-party entanglement and demonstrate the existence of a function whose computation requires strictly lesser classical communication complexity compared to the scenario where no quantum entanglement is used.

Singh, Kumar and Pal  show how $`n`$ nodes connected in a network can be entangled together in a maximally entangled $`n`$-partite state $`(|0_10_2\mathrm{}0_n>+|1_11_2\mathrm{}1_n>)\sqrt{2}`$ using linear classical communication complexity provided the nodes were earlier entangled in $`(n1)`$ EPR pairs. The $`(n1)`$ EPR pairs given initially are required to be shared between agents forming edges of a spanning tree in the network; such a spanning tree is called the spanning EPR tree . Brassard et al. show in  that prior multi-partite entanglement can be used by $`n`$ agents in a perfect quantum protocol (not using any form of communcation between the $`n`$ agents), to solve a multi-party distributed problem, whereas, no classical deterministic protocol succeeds in solving the same problem with probability away from $`\frac{1}{2}`$ by a fraction that is larger than an inverse exponential in the number of agents. Multi-partite entanglement is also used by Buhrman, van Dam, Hoyer and Tapp  where an $`n`$-party problem was shown to have quantum communication complexity $`n`$ qubits, whereas the classical communication complexity was shown to be $`\mathrm{\Omega }(n\mathrm{log}n)`$ cbits.

In this paper we show that a sufficient amount of multi-partite quantum entanglement helps in electing a leader unambiguously in a distributed network of agents where it is not posible for any agent to introduce any statistical bias in the election. We show in our Protocol I that a total of $`O(\mathrm{log}n)`$ distinct multi-partite entanglement states are capable of supporting a protocol with linear classical communication complexity for leader election even in the presence of nodes that may lie or be biased. Here, $`n`$ is the number of nodes in the network. Unbiased election is enabled by the sharing of multiple multi-partite entangled states. Using one honest referee, we also ensure that all honest nodes agree on the elected leader. We also demonstrate the difficulty of performing unbiased and fair election of a leader with linear classical communication complexity in the absence of quantum entanglement, even if all nodes have perfect random bit generators. We demonstrate this through our Protocol II where we show the following separation result: too many biased liars would prevent the referee from conducting a fair and unbiased election, whereas, a small number of liars would matter little. This separation result has two aspects. Firstly, it shows as expected that the probability of the election being biased is small and indeed vanishes as the number of agents grows asymptotically when there is a small number of biased agents. So, we observe that it is possible to conduct fair and unbiased election in our randomized Protocol II with small error probability with linear communication cost. Protocol I performs fair and unbiased election but needs $`O(\mathrm{log}n)`$ multi-partite maximally entangled states or $`ebits`$. Note however that even though the number of biased voters ($`\mathrm{log}n`$) is small, it is an increasing function of $`n`$. So, with small but growing number of liars, the randomized Protocol II achieves vanishingly small biased election probability as $`n`$ grows asymptotically; this demonstrates partial tractability of our leader election problem with pure randomization and no quantum resource. Secondly, the finite limiting probability $`e^C,C>0`$ ($`C`$ is a constant), of fair and unbaised election with the number of biased voters as high as $`n/\mathrm{log}n`$ (Protocol II) shows that randomization cannot handle too many cheating agents who are essentially biased in a statistical sense. Even here, we keep the large number of biased agents $`o(n)`$ by choosing $`n/\mathrm{log}n`$ biased agents.

Further, we define two complexity classes (i) $`QECC(g,n)`$, the class of problems solvable using $`g`$ $`ebits`$ shared amongst $`n`$ agents with $`O(n)`$ cbit communication complexity, and, (ii) $`RCC(r,n)`$, the class of problems solvable using $`r`$ random bits shared bewteen $`n`$ agents, with $`O(n)`$ cbit communication complexity. We show that $`RCC(\mathrm{log}n,n)QECC(\mathrm{log}n,n)`$. The problem of determining whether this is a proper containment is an interesting fundamental open problem concerning the power of quantum entanglement as a computational and communication resource.

Throughout our analysis, we follow a uniform statistical notion of fairness and bias, as we define shortly below. In short, an unbiased agent is supposed to always vote a 0 or a 1 with equal probability as per an ideal random bit generator that the agent possesses. So, any biased agent, using differing probabililties for 0 and 1 votes, would cause some bias towards some agent addresses in the election of the leader. Fairness is ensured by giving agents equal opportunity in voting for determination of the leader.

On the one hand, the successful completion of a protocol would require all agents in the network to agree on an elected leader. On the other hand, agents may be biased in voting for a leader. It is therefore also a matter of concern that one agent amongst all in the network must act as a referee, trying to conduct a fair and unbiased election process. Since election essentially requires voting by all agents in a fair setup and, agreement by all agents about the elected leader, we observe that linear classical communication complexity is essentially a lower bound in our model. In this sense, our model differs from the models of Yao  and also those of .

In Section 2 we introduce our own definitions and notions of fairness, bias, referee and agreement. We essentially propose two leader election protocols, both use only linear amounts of classical communication. The first one (Protocol I in Section 3) assumes the presence of a sufficient number of $`n`$-partite maximally entangled states, and, the second one (Protocol II in Section 4) assumes no such entanglement but only perfect random bit generation capability in each agent. We show a separation result by demonstrating that too many biased liars would prevent the referee from conducting a fair and unbiased election, whereas, a small number of liars would matter little. In Section 5 we define two new complexity classes related to the problem of unbiased election of a leader and pose a few open problems.

## 2 The Communication Framework

We wish to have a complete protocol with a starting signal going to each agent, signifying initiation. This requires communicating a constant number of bits of classical information (cbits) from one special agent to the rest. (There may be a constant number of special tasks like the starting and ending of protocols. So a constant number of bits are sufficient for such signals.) We fix one arbitrary node to initiate the protocol by sending this initial message to all other agents. This node is also set as the referee node; the referee is supposed to be conducting the protocol and collecting agreements from all agents, including itself, about the elected leader. If an agent does not accept the elected leader, he is simply left out; further computation can be carried out by only those agents who accept the election of the leader. In other words, we require total agreement over the elected leader by all the agents if all of them wish to continue computations together, after leader election. Such an agreement comes in the form of a cbit from each agent after the leader is elected. The referee first informs the leader about his election and the leader then sends a constant number of cbits to each agent, seeking agreement or approval.

Identifying a leader amounts to fixing $`\mathrm{log}n`$111All logarithms in this paper are to the base 2. bits to uniquely address an agent out of $`n`$ agents. Any protocol determining a leader must therefore come up with these bits. We assume without loss of generality that $`n`$ is a power of 2 and therefore we need to determine exactly $`\mathrm{log}n`$ bits.

Now we state the communication and message passing model for classical communication in our protocols. We assume that units of information (cbits) are communicated from a source agent $`S`$ to a destination agent $`T`$. The destination can store as many addresses of senders of incoming messages so that in future it may choose to reply to any of them. We also assume that the cost of communicating a single bit of information between any source $`S`$ and any destination $`T`$ is the same and is counted as one cbit, independent of the network connecting the agents. We are now in a position to derive our protocols.

We define bias as follows. If an agent $`A_i`$ can somehow influence the election of the leader in a way that is statistically different from the way another agent $`A_j`$ influences election, then we say that there is a bias in the election process. In other words, the determination of the address of the leader must be equally influenced by each participating agent.

The notion of fairness is subtlely distinct from that of bias. Here, fairness means ensuring that all agents get a fair deal in participating in the election process, in a statistical sense. These notions of statistical bias and fairness will become more explicit when we analyse our protocols in the next two sections.

## 3 Fair and Unbiased Leader Election Using $`n`$-partite Entanglement

We assume that the $`n`$ agents $`A_1,A_2,\mathrm{},A_n`$ are in multiple $`n`$-partite maximally entangled states presented in the computational bases. Such states are called $`ebits`$ . Each agent houses his own qubits in these entangled states. We also assume that there is a sufficiently large number of such states. In particular, we choose to have $`\mathrm{log}n`$ such $`ebit`$ states $`S_1,\mathrm{},S_{\mathrm{log}n}`$, with each agent $`A_i,1in`$ housing $`\mathrm{log}n`$ qubits $`qu_{i,j},1j\mathrm{log}n`$, one for each of the $`\mathrm{log}n`$ states. The state $`S_k`$ is therefore $`(|0_{1,k}0_{2,k}\mathrm{}0_{n,k}>+|1_{1,k}1_{2,k}\mathrm{}1_{n,k}>)/\sqrt{2}`$ for $`1k\mathrm{log}n`$, where the subscript $`k,i`$ denotes the $`k`$-th qubit of the $`i`$-th agent. Clearly, a total of $`n\mathrm{log}n`$ qubits are used. It is possible to create such a maximally entangled multi-partite (ebit) state using a protocol in  that requires linear classical communication complexity and uses up $`(n1)`$ EPR pairs generated apriori in the form of an EPR spanning tree of the network. Our Protocol I is initiated by the referee who sends one cbit to all agents signalling the start of the protocol. On getting this trigger, each node $`A_i,1in`$ performs measurements on its $`\mathrm{log}n`$ qubits $`qu_{i,j},1j\mathrm{log}n`$. These measurements are performed in the respective computational bases for each of the $`\mathrm{log}n`$ entanglement sets. Whoever measures a qubit, gets a value zero or one, with probability $`\frac{1}{2}`$ for each of 0 and 1 values. Indeed, this is a perfect random bit generator. Also, no matter in what order $`A_i`$ and $`A_j`$, $`1j<in`$ measure their $`k`$-th qubits, both of them get the same answer. This follows from the postulates of quantum mechanics and the property of maximal entanglement for state $`S_k`$. The referee can now use these $`\mathrm{log}n`$ measured bits to produce an address or index $`m`$ for the agent $`A_m,1mn`$ and accept it as an unbiased leader. Note that such a leader is perfectly random and therefore totally unbiased. There is no hiding of the measurements made by the agents about their $`k`$-th qubits since the referee can always check the measured value of the $`k`$-th qubit of the $`k`$-th entangled state $`S_k`$ by reading its own $`k`$-th qubit. Even if the referee seeks this value from another agent and the agent lies, the referee can simply ignore any false information provided because it knows the true value of the measured $`k`$-th qubit. Now it is the turn of all agents to send one cbit to the referee and express their agreement over the elected leader. If an agent agrees on the leader by sending a cbit 1, it is fine and the agent henceforth sticks to the commitment of having agreed on the leader. Otherwise, we say that the agent is discarded and only agreeing agents are considered to contuinue meaningful computations any further. We also assume that the elected leader is committed to act as leader henceforth. We summarize the lack of bias and the guarantee of fairness in the election of the leader in the above Protocol I in the following theorem.

###### Theorem 3.1

It is possible to conduct an unbiased and fair election of a leader amongst $`n`$ agents who share $`\mathrm{log}n`$ maximally entangled $`n`$-partite states with $`O(n)`$ bits of classical communication.

## 4 Randomized Election

We believe that randomization alone is not capable of supporting distributed protocols that can achieve fair and unbiased leader election provided only a linear number of cbits are exchanged in the protocol. Although multi-partite entanglement helps us perform fair and unbaised leader election with only linear classical communication complexity (as we see in Protocol I in the previous section), randomization alone, coupled with only linear number of cbits, may lead to fair and unbiased election of a leader only in certain special and restricted cases. Recall that the leader selected in the case of Protocol I is perefctly random, irrespective of the order in which the agents make measurements on the entangled qubits; also, each agent is equally likely to become the leader. In this section we attempt to mimic similar unbiased and fair election and we propose a Protocol II which succeeds in certain cases in a randomized sense. To this effect, we design our Protocol II by making the agents together play the role of an adversary to the referee; the leader declared by the referee may be biased if an agent so chooses. The evil designs of the agents may however not work always and therefore a finite probability of fair and unbiased election may result.

What really is fairness and bias in this scenario, respecting the overall definition we gave earlier in Section 2? Firstly, fairness forces us to design protocols that treat all agents equally, providing them with equal opportunity in the election process, statistically speaking. This must be ensured at the referee’s end. Secondly, bias can be introduced in the choice presented by an agent so that individual agents have unequal probabilities of becoming elected as leader. So, for unbiased election, we may expect agents to honestly use their ideal random bit generating capability. A biased agent (adversary) may not use its perfect random bit generator and may instead bias all his outcomes in a protocol. Providing such biased bits to the referee would lead to a biased election where all the agents may not have equal probability of getting elected. So, we propose our Protocol II where each agent sends just one bit (biased or perfectly ideal random bit) to the referee. Being honest and fair, the referee collects all these $`n`$ bits from the $`n`$ agents and then selects a perfectly random sample of $`\mathrm{log}n`$ bits from the $`n`$ collected bits to compute the index $`m`$ of the leader $`A_m`$. Initiation of the protocol and the agreement step are as in Protocol I. This completes the description of Protocol II.

We now show that biased adversary agents can not fool the referee in an asymptotic and statistical sense, if there are $`k=\mathrm{log}n`$ biased agents. We also show that a larger number of biased agents, proportional to $`n/\mathrm{log}n`$, can turn the table the other way; in this case the probability that the referee can not vanquish the evil designs of the biased agents is a nonzero constant. We state this as a theorem.

###### Theorem 4.1

It is possible to design a randomized protocol that uses no quantum entanglement, but uses perfect random bit generating capability in each of the $`n`$ agents and only $`O(n)`$ bits of classical communication so that as $`n`$ grows asymptotically, (i) fair and unbiased election of a leader is guaranteed with probability approaching 1 if there are only $`O(\mathrm{log}n)`$ biased agents, and (ii) unbiased election is possible only with a nonzero finite limiting probability if the number of biased agents is $`O(n/\mathrm{log}n)`$.

Proof. The referee collects the single bit votes sent by the agents, each agent sending one bit vote. To be fair to all, the sincere referee picks up a random sample of $`\mathrm{log}n`$ bits from the $`n`$ votes and arranges the bits in random order to create the $`\mathrm{log}n`$ bit address of the leader who is thus deemed elected. Note here that as long as each agent is honestly sending a vote which is a perfect random bit, the address generated for the leader will be fair and statistically uniform over all agents. So, as long as the agents themselves are honest, leader election is fair even about who becomes the leader. The problem arises when the bit sent by a dishonest or biased agent is not perfectly random. If the referee happens to select any such biased bit in his random sample for the leader’s address, then the leader elected is statistically biased, not uniform over all agents. We note that the probability that the random sample would be able to avoid biased bits is $`p(n,k)=C(nk,\mathrm{log}n)/C(n,\mathrm{log}n)`$ if there are $`k<n\mathrm{log}n`$ biased bits. We show in Appendix that $`p(n,\mathrm{log}n)`$ approaches 1 as $`n`$ grows asymptotically. Thus, we conclude that as long as there are only a small number of liars amongst the agents who attempt to bias leader election by casting biased votes, the referee succeeds in deciding upon a leader without bias with very high probability approaching unity. We also show in Appendix that $`p(n,n/\mathrm{log}n)`$ approaches a non-zero constant as $`n`$ grows asymptotically, showing that there is a non-zero probability that the referee fails to deliver unbiased election when the number of liars biasing their votes is as large as $`O(n/\mathrm{log}n)`$.

## 5 New Complexity Classes

In the context of our results (Theorems 3.1 and 4.1), we define two new complexity classes where the $`n`$ agents are required to use at most linear classical communication but may use some quantum resource or randomized resource. We say that a problem belongs the class $`QECC(g,n)`$ if the agents use $`O(g)`$ $`ebits`$ and $`O(n)`$ $`cbits`$. An $`ebit`$ is equivalent to an $`n`$-partite maximally entangled state shared between the agents as we have used in Section 3. We saw in our Theorem 3.1 how the fair and unbiased election of a leader is possible in $`QECC(\mathrm{log}n,n)`$. We define another class $`RCC(r,n)`$, of problems that can be solved using $`n`$ agents that share an $`O(r)`$-bit random string and use only linear classical communication complexity. It is easy to observe that $`RCC(r,n)QECC(r,n)`$. This is because by measuring his own qubits, one qubit in each of the ebits shared with other agents, an agent can generate as many ideal random bits shared between all the agents. Whether this relationship is a proper containment is an interesting fundamental open problem concerning the power of quantum entanglement as a computational and communication resource.

## 6 Concluding Remarks

It would be interesting to design other classical and randomized protocols for fair and unbiased leader election that work for any number of liars, possibly using superlinear communication complexity. Other fundamental problems in distributed computing related to multicasting, broadcasting or consensus formation may also be studied to see if quantum entanglement can reduce classical communication requirement in asymptotically randomized protocols.

## Appendix

Let $`C(n,r)`$ be the number of combinations of $`n`$ things taken $`r`$ at a time without repetitions. We note that

$$\frac{C(nk,\mathrm{log}n)}{C(n,\mathrm{log}n)}=\frac{(nk)!/(\mathrm{log}n)!(nk\mathrm{log}n)!}{n!/(\mathrm{log}n)!(n\mathrm{log}n)!}$$

Cancelling the common terms and applying Sterling’s approximation to the four terms for large $`n>>k`$, we have the estimate

$$\frac{(nk)^{nk+\frac{1}{2}}(n\mathrm{log}n)^{n\mathrm{log}n+\frac{1}{2}}}{(nk\mathrm{log}n)^{nk\mathrm{log}n+\frac{1}{2}}n^{n+\frac{1}{2}}}$$

$$=\frac{n^{nk+\frac{1}{2}}(1\frac{k}{n})^{nk+\frac{1}{2}}(n\mathrm{log}n)^{n\mathrm{log}n+\frac{1}{2}}}{(n\mathrm{log}n)^{nk\mathrm{log}n+\frac{1}{2}}(1\frac{k}{n\mathrm{log}n})^{nk\mathrm{log}n+\frac{1}{2}}n^{n+\frac{1}{2}}}$$

$$=\frac{(1\frac{\mathrm{log}n}{n})^k(1\frac{k}{n})^{nk+\frac{1}{2}}}{(1\frac{k}{n\mathrm{log}n})^{nk\mathrm{log}n+\frac{1}{2}}}$$

Now note that putting $`k=\mathrm{log}n`$, and using Euler’s formula with $`n`$ increasing asymptotically, the above estimate approaches unity. On the other hand, putting $`k=\frac{Cn}{\mathrm{log}n}`$, for a constant $`C>0`$ takes the estimate to $`e^C`$ as $`n`$ blows up. Here $`e`$ is Euler’s constant.


text:
# First normal stress difference and crystallization in a dense sheared granular fluid1footnote 11footnote 1 Physics of Fluids, vol. 15, no. 8 (2003)

## I Introduction

In the last decade, a lot of research activity took place to unveil the properties of granular materials1,2, primarily because of their industrial importance, but also due to their fascinating properties. This has unraveled many interesting and so far unresolved phenomena (for example, clustering, size-segregation, avalanches, the coexistence of gas, liquid and solid, etc.). Under highly excited conditions, granular materials behave as a fluid, with prominent non-Newtonian properties, like the normal stress differences3. While the normal stress differences are of infinitesimal magnitudes in a simple fluid (e.g. air and water), they can be of the order of its isotropic pressure in a dilute granular gas4. From the modelling viewpoint, the presence of such large normal-stress differences readily calls for higher-order constitutive models5,6 even at the minimal level.

Studying the non-Newtonian behaviour is itself an important issue, since the normal stresses are known to be the progenitors of many interesting and unique flow-features (e.g. rod-climbing or Weissenberg-effect, die-swelling, secondary flows, etc.7) in non-Newtonian fluids. Also, normal stresses can support additional instability modes (for example, in polymeric fluids and suspensions7-10, which might, in turn, explain some flow-features of granular fluids. For example, particle-clustering11-13 has recently been explained from the instability-viewpoint using the standard Newtonian model for the stress tensor12,14,15

The kinetic theory of Jenkins & Richman16 first showed that the anisotropy in the second moment of the fluctuation velocities, due to the inelasticity of particle collisions, is responsible for such normal stress behaviour. They predicted that the first normal stress difference (defined as $`𝒩_1=(\mathrm{\Pi }_{xx}\mathrm{\Pi }_{yy})/p`$, where $`\mathrm{\Pi }_{xx}`$ and $`\mathrm{\Pi }_{yy}`$ are the streamwise and the transverse components of the stress deviator, respectively, and $`p`$ is the isotropic pressure, see section IIB) is maximum in the dilute limit, decreases in magnitude with density, and eventually approaches zero in the dense limit. Goldhirsch & Sela4 later showed that the normal stress differences appear only at the Burnett-order-description of the Chapman-Enskog expansion of the Boltzmann equation. Their work has clearly established that the origin of this effect (in the dilute limit) is universal in both atomic and granular fluids, with inelasticity playing the role of a magnifier and thus making it a sizeable effect in granular fluids. While the source of the normal stress differences in the dilute limit has been elucidated both theoretically and by simulation, its dense counterpart has not received similar attention so far. This is an important limit since the onset of dilatancy (volume expansion due to shear17,18), crystallization, etc. occur in the dense regime, which in turn would influence the normal stress differences.

Previous hard-sphere simulations19,3,11 did look at the normal stress differences, but they did not probe the dense limit in a systematic way. These simulations showed that the first normal stress difference vanishes in the dense limit, in line with the theoretical predictions of Jenkins & Richman16. On the other hand, the soft-sphere simulations of Walton & Braun20, with frictional particles, showed that this quantity can change sign in the same limit. Our work with smooth inelastic hard-disks unequivocally demonstrates that $`𝒩_1`$, indeed, changes its sign at some critical density in the dense regime, due to the sign-change of its collisional component at a critical density, which depends on the value of the coefficient of restitution ($`e`$). More importantly, we show that the origin of $`𝒩_1`$ in the dense limit is distinctly different from that in a dilute granular gas. At the microstructural-level, certain topological changes in the anisotropic structure of the collision-angle distribution with density are responsible for the observed sign-reversal of $`𝒩_1`$.

We use the familiar smooth hard-disk model for an event-driven simulation21 of the uniform shear flow configuration, focussing mainly on the normal stress behaviour and the microstructure formation as functions of the density and inelasticity. The details of the simulation technique and the relevant macroscopic quantities are described in section II. The simulation results on the first normal stress difference, the radial distribution function, the collision angle distribution and the crystalline-structure are presented in section III. Possible modelling approaches to incorporate the normal stress differences are discussed in section IV. In section V we summarize our findings, with suggestions for possible future work.

## II Simulation method

We consider a collection of smooth inelastic hard-disks in a square box of size $`\stackrel{~}{L}`$ under uniform shear flow — let $`\stackrel{~}{x}`$ and $`\stackrel{~}{y}`$ be the streamwise and transverse directions, respectively, with the origin of the coordinate-frame being positioned at the centre of the box. The snapshot of a typical simulation, with non-dimensional coordinates, is shown in Fig. 1($`a`$). Note that the dimensional quantities are denoted by tildes, and the reference length, time and velocity scales for non-dimensionalization will be specified later in this section.

Let the diameter and the mass of the particle be $`\stackrel{~}{\sigma }`$ and $`\stackrel{~}{m}`$, respectively. The pre- and post-collisional particle velocities of particle $`1`$ are denoted by $`\stackrel{~}{𝐜}_1`$ and $`\stackrel{~}{𝐜}_1^{}`$, respectively. Hence, the velocity of particle $`2`$ relative to $`1`$ is $`\stackrel{~}{𝐜}_{21}=\stackrel{~}{𝐜}_2\stackrel{~}{𝐜}_1`$. Let $`𝐤_{21}=𝐤`$ be the unit vector directed from the center of particle $`2`$ to that of particle $`1`$ at contact. The pre- and post-collisional velocities are related by the expression:

$$𝐤\stackrel{~}{𝐜}_{21}^{}=e(𝐤\stackrel{~}{𝐜}_{21}),$$
(1)
where $`e`$ is the coefficient of normal restitution, with $`0e1`$; note that we restrict ourselves to perfectly smooth particles. The expression for the collisional impulse is

$$\stackrel{~}{𝐈}=\stackrel{~}{m}(\stackrel{~}{𝐜}_1^{}\stackrel{~}{𝐜}_1)=\frac{\stackrel{~}{m}}{2}(1+e)(𝐤\stackrel{~}{𝐜}_{21})𝐤,$$
(2)
directed along $`𝐤`$.

### II.1 Model system and algorithm

The system is periodic in $`\stackrel{~}{x}`$-direction, i.e. a particle crossing the left/right boundary re-enters the system through the opposite boundary at the same vertical position with unchanged velocities. To impose a uniform shear rate ($`\stackrel{~}{\gamma }=\stackrel{~}{U}/\stackrel{~}{L}`$) in the $`\stackrel{~}{y}`$-direction, the top and bottom image boxes, bounding the central box, are set in motion with velocities $`\stackrel{~}{U}/2`$ and $`\stackrel{~}{U}/2`$, respectively, in the streamwise direction. This is the standard approach to attain the state of uniform shear flow (USF) by imparting momemtum transfer by shearing, originally introduced by Lees & Edwards22. Overall, this system represents an extended doubly-periodic system where the periodicity in the transverse direction is in the local Lagrangian frame. In a typical simulation, the disks are initially placed randomly in the computational box, and the initial velocity field is composed of the uniform shear and a small Gaussian random part. An event-driven algorithm is then used to update the system in time, the details of which may be found in Alam & Luding23,24.

To ascertain whether the system has reached the statistical steady-state, the time evolution of the average fluctuation kinetic energy (‘granular’ energy, defined in the next section) is monitored, see Fig. 1($`b`$). Due to the balance between the shear work and the collisional dissipation under homogeneous shear deformation, the granular energy attains a constant value at the steady state. Depending on the value of the coefficient of restitution and the number of particles, it takes about thousand collisions per particle to reach such a statistical steady-state – the lower the value of $`e`$, the more quickly the system reaches the steady state and vice versa. The simulation is then allowed to run for another $`15000`$ collisions per particle to gather data to calculate the macroscopic quantities. A few longer runs (30000 collisions per particle) were also checked, with no significant change on the measured quantities. Another quantity which was simultaneously monitored, along with granular energy, is the linearity of the streamwise velocity profile across the Couette gap, and we found that the calculated shear rate (i.e. the slope of the velocity profile) fluctuated around the imposed shear rate by at most $`1\%`$ at densities where crystallization is not evidenced.

### II.2 Macroscopic quantities

With $`\stackrel{~}{L}`$, $`\stackrel{~}{\gamma }^1`$, $`\stackrel{~}{\gamma }\stackrel{~}{L}`$, and $`\stackrel{~}{m}`$, used as the reference length, time, velocity, and mass, respectively, the relevant dimensionless quantities are:

$$\sigma =\frac{\stackrel{~}{\sigma }}{\stackrel{~}{L}},(𝐜,𝐮,𝐂)=\frac{1}{\stackrel{~}{\gamma }\stackrel{~}{L}}(\stackrel{~}{𝐜},\stackrel{~}{𝐮},\stackrel{~}{𝐂}),𝐏=\frac{\stackrel{~}{𝐏}}{\stackrel{~}{\rho }\stackrel{~}{\sigma }^2\stackrel{~}{\gamma }^2},T=\frac{\stackrel{~}{T}}{\stackrel{~}{\sigma }^2\stackrel{~}{\gamma }^2},$$
(3)
where $`𝐮`$ is the ‘hydrodynamic’ velocity, $`𝐂=𝐜𝐮`$ the fluctuation (peculiar) velocity of particles, $`\stackrel{~}{\rho }`$ the material density of particles, $`𝐏`$ the stress tensor, and $`T`$ the granular energy.

The macroscopic stress, as measured in discrete particle simulations, is a byproduct of the particle-level mechanisms of momemtum transfer. As in the hard-core model of dense gases, the stress is the sum of its kinetic and collisional components. The former arises from the transport of momentum as the particles move through the system carrying their momentum, while the latter is due to the direct interparticle collisions. The homogeneity of the uniform shear flow allows us to calculate the stress by averaging it over the whole computational box3,23,24.

The stress, defined in the compressive sense, may be decomposed in the standard way:

$$𝐏=𝐏^k+𝐏^c=p\mathrm{𝟏}+𝚷,$$
(4)
where $`p`$ is the pressure, $`𝚷`$ the pressure deviator and $`\mathrm{𝟏}`$ the unit tensor. From the off-diagonal components of the pressure deviator, we can calculate the shear viscosity which relates the rate of strain to the shear stress:

$$\mu =\mathrm{\Pi }_{xy}/\frac{\mathrm{d}u}{\mathrm{d}y}.$$
(5)
For the steady uniform shear flow, thus, the dimensionless shear viscosity can also be interpreted as the shear stress due to our adopted scaling, $`\mathrm{d}u/\mathrm{d}y=\gamma =1`$. The diagonal components of the pressure deviator can be non-zero, giving rise to normal stress differences. The first normal stress difference is defined as

$$𝒩_1=\frac{(\mathrm{\Pi }_{xx}\mathrm{\Pi }_{yy})}{p}.$$
(6)
Note that we have scaled this quantity by pressure to ascertain its relative magnitude with respect to pressure. For a standard Newtonian fluid, $`𝒩_1=0`$ and thus $`𝒩_1`$ is an indicator of the non-Newtonian character of the fluid. $`𝒩_1`$ can be decomposed into kinetic and collisional parts:

$$𝒩_1=𝒩_1^k+𝒩_1^c=\frac{(\mathrm{\Pi }_{xx}^k\mathrm{\Pi }_{yy}^k)}{p}+\frac{(\mathrm{\Pi }_{xx}^c\mathrm{\Pi }_{yy}^c)}{p}.$$
(7)
Note that the sign of $`𝒩_1`$ crucially depends on the convention used to define the stress tensor. For example, in the rheology literature, stress is typically defined in the tensile sense7. A positive $`𝒩_1`$ for the compressive case is equivalent to its negative value for the tensile case and vice versa. This point should be kept in mind while making any comparison with data in the rheology literature.

From the trace of the kinetic stress tensor, $`𝐏^k`$, one can calculate the granular energy,

$$T=\frac{1}{2\sigma ^2}\left[\frac{1}{N}\underset{i=1}{\overset{N}{}}C_iC_i\right],$$
(8)
which is a measure of the random motion of the particles with respect to the mean motion.

There are two dimensionless control parameters: the volume fraction of particles ($`\nu `$) and the coefficient of normal restitution ($`e`$). The shear rate is also a control parameter, however, due to normalization we have $`\gamma =\stackrel{~}{\gamma }/\stackrel{~}{\gamma }=1`$, and changing the value of $`\stackrel{~}{\gamma }`$ does not influence the reported results; the imposed shear rate $`\gamma `$ is thus kept fixed at unity. The simulations are carried out for the whole range of solid volume fractions, varying from the dilute to the dense limit, over a large range of values for the coefficient of restitution ($`e=0.3`$$`0.99`$). For most of the simulations, the number of particles are fixed to $`N=1024`$, and increasing the value of $`N`$ by fourfold ($`N=4096`$) did not affect the reported quantities noticeably; for example, the change in $`𝒩_1`$ was about $`4.3\%`$ and $`5.1\%`$ at $`\nu =0.3`$ and $`0.75`$, respectively, for a restitution coefficient of $`e=0.9`$. We note here that the system-size dependence of the rheological quantities (pressure and viscosity) is known to be strong only for a small number of particles ($`N<100`$)11,24,25.

For the typical simulation in Fig. 1($`a`$), at steady-state, after $`2\times 10^7`$ collisions, the parameter values were $`\nu =0.5`$, $`N=1024`$ and $`e=0.7`$. The variations of the granular energy $`T`$ and the calculated shear rate $`\gamma _{cal}`$ with time are shown in Fig. 1($`b`$), along with corresponding initial variations in two insets. (The data represent the instantaneous values of $`T`$ and $`\gamma _{cal}`$ sampled at a regular interval of $`400`$ collisions – no time averaging is involved here.) Note that $`\gamma _{cal}`$ was computed by binning the system into $`20`$ equal-size bins in the transverse direction and taking averages over all particles in each bin. It is observed that the granular energy reaches its steady value ($`T=0.6121\pm 0.022`$) quickly after the initial transients and the calculated shear rate fluctuates around its imposed value ($`\gamma =1`$) by about $`\pm 1\%`$. The fluctuations in both $`T`$ and $`\gamma _{cal}`$ at steady-state are due to the finite-size of the system and diminish with increasing number of particles as $`N^{1/2}`$.

## III Results

For detailed results on the quantities pressure, shear viscosity and granular energy, and for their comparison with kinetic theory predictions, we refer to our recent study24. Here, we will mainly focus on the behaviour of the first normal stress difference and its kinetic and collisional components. We also present results on the pair distribution function and the collision angle distribution to characterize microstructures. Lastly we will present results on crystal-formation at high densities, signalling the coexistence of fluid and solid, complementing recent results in non-sheared systems26,27.

### III.1 Normal stress difference

Figure 2($`a`$) shows the variation of the first normal stress difference ($`𝒩_1`$) with density for two values of the coefficient of restitution. It is observed that $`𝒩_1`$ is maximum at the dilute limit and decreases thereafter with $`\nu `$. The overall variation of $`𝒩_1`$ with density looks similar at other values of the coefficient of restitution, with a difference in the magnitude of $`𝒩_1`$. The inset in Fig. 2($`a`$) shows that $`𝒩_1`$ decreases quite sharply in the dense limit and becomes negative at some density ($`\nu =\overline{\nu }`$). Increasing the value of the restitution coefficient decreases this critical density $`\overline{\nu }`$. The arrows on the left-ordinate indicate the asymptotic values of $`𝒩_1`$ for a two-dimensional granular gas in the dilute limit4:

$$𝒩_1=1.0448(1e^2).$$
(9)
The anisotropy in the second moment of the fluctuation velocity is primarily responsible for the finite normal stress difference in the dilute limit16,28 and this shows up only at the Burnett-order of the Chapman-Enskog expansion4. We should mention here that the limit $`e1`$ is singular and the normal stress difference survives even in the elastic limit as pointed out by Goldhirsch & Sela4. The corresponding expression for $`𝒩_1`$ in a molecular gas is:

$$𝒩_11.358\frac{\gamma ^2\mathrm{}^2}{u^2},$$

where $`\mathrm{}`$ is the mean free path and $`u^2`$ is the rms of the velocity fluctuations. However, because of its extremely small magnitude under normal conditions, the normal stress difference is not measurable in a molecular fluid.

Previous hard-sphere simulations of Campbell and coworkers3,19 are in variance with our result in that they found $`𝒩_10`$ as $`\nu \nu _{max}`$. However, the soft-sphere simulations of Walton & Braun20, with frictional particles, support our observation that $`𝒩_1`$ indeed undergoes a sign-reversal. To better understand what is responsible for the sign-reversal of $`𝒩_1`$, we look at the kinetic and collisional components of the first normal stress difference. Figure 2($`b`$) shows the variations of $`𝒩_1^k`$ and $`𝒩_1^c`$ with density at $`e=0.7`$. We observe that $`𝒩_1^k`$ is maximum at the dilute limit and decreases monotonically to zero as $`\nu `$ approaches the packing limit. Except for the dense limit, the overall behaviour of $`𝒩_1^k`$ represents that of the total normal stress difference. The collisional component, $`𝒩_1^c`$, shows a non-monotonic variation with density: $`𝒩_1^c`$ is zero in the dilute limit, increases with increasing $`\nu `$, remains almost constant for intermediate densities, and then decays sharply in the dense limit. Interestingly, $`𝒩_1^c`$ becomes negative at some critical density ($`\nu =\overline{\nu }`$) beyond which the behaviour of $`𝒩_1^c`$ mirrors that of $`𝒩_1`$ (see inset). Thus the normal stress behaviour in the dense regime is clearly due to the anisotropy in the collisional stress.

Recall that the kinetic theory of Jenkins & Richman16,28 predicts that $`𝒩_10`$ in the dense limit. The predictions of the revised Enskog theory of Santos et al.29 are in line with that of Jenkins & Richman, even though their kinetic model is claimed to be valid even in the crystalline-phase. Since the source of normal stress differences in all these theories is linked to the anisotropy in the second moment of velocity fluctuations (which vanishes as $`\nu \nu _{max}`$), they are unable to predict the correct behaviour of normal stresses in the dense limit. We would like to stress here that, as mentioned in the Introduction, many fascinating non-Newtonian effects7 are primarily determined by the first normal stress difference and its sign. For example, the rod-climbing effect can occur in a granular fluid if $`𝒩_1<0`$, i.e. only in the dense limit (assuming that the second normal stress difference is negligible). Furthermore, having a constitutive model which reproduces the correct sign of $`𝒩_1`$ is also important, since it is well known in the rheology literature that the odd sign for $`𝒩_1`$ leads to the instability of the rest state30.

The sign-reversal of $`𝒩_1^c`$ can be succinctly presented as a phase-diagram in the ($`\nu ,e`$)-plane by plotting the zeros of $`𝒩_1^c`$ as a function of the coefficient of restitution, see Fig. 3. Below the solid line, $`𝒩_1^c`$ is positive, and negative above it. Also plotted in this figure is the line for the zeros of $`𝒩_1`$ which, as expected, lies slightly above. Thus, the first normal stress difference is zero along the solid line which may be called the symmetry-line. It is observed that decreasing the coefficient of restitution increases the critical density ($`\overline{\nu }`$) at which $`𝒩_1^c`$ changes sign. As we approach the elastic limit, $`\overline{\nu }`$ depends strongly on the value of $`e`$. We further note that as $`e1`$, $`\overline{\nu }0.62`$ which is well below the freezing-point density of a 2D hard-disk fluid, $`\nu _f0.70`$26,27.

We need to mention here that the sign-change of $`𝒩_1`$ is not uncommon in other non-Newtonian fluids. For example, in non-Brownian viscous suspensions, $`𝒩_1`$ changes sign at high Peclet number31,32. However, the reason for this effect is quite different in granular fluids as we show below.

### III.2 Microstructural features

To understand the microstructural mechanism for the origin of the first normal stress difference and its sign-reversal, here we probe several microstructural features of a dense granular fluid.

#### III.2.1 Radial distribution function

Typical snapshots of the system in the dense regime are shown in Fig. 4 at four different densities with $`e=0.7`$. Note that the density for the subplot 4($`c`$) is $`\nu =0.725`$ for which $`𝒩_1^c0`$. Looking at the corresponding distribution of granular energies (not shown here for brevity), we could find signatures of clusters (group of particles) with lower energies surrounded by particles with higher energy. To understand the flow-microstructures and their energetics at such high densities, we need to probe the pair distribution function and similar measures for the structure of the packing.

Figure 5($`a`$-$`d`$) shows the radial distribution function $`g(r)`$ at four densities, with parameter values as in Fig. 4. The thin, dotted lines are data from a non-sheared, homogeneous, elastic system26,27, whereas the thick, solid lines represent a sheared situation with rather strong dissipation $`e=0.7`$. (Note that the elastic distribution function has been measured under non-sheared periodic boundary conditions; also for $`e`$ somewhat smaller than unity, the same results were obtained as long as the system remains homogeneous.) We observe that the weak difference at low density $`\nu =0.6`$ grows with increasing density, concerning two aspects:
(i) The peak value of contact in the sheared systems is always larger than that in a homogeneous system of the same density, and the difference increases strongly with density, another indicator for clustering33.
(ii) The peaks and valleys, which allow to distinguish between different lattice structures, are different in the sheared case when compared to the non-sheared situation. In the former case, peaks at $`r/\sigma =1,2,3,\mathrm{}`$ are observed, signalling shell-formation about any test-particle. In the latter case, the peaks at $`r/\sigma =1,\sqrt{3},2,\mathrm{}`$ indicate a crystallization transition and the development of a triangular lattice.

The peaks in the sheared situation, e.g. at $`r/\sigma =2,3`$, become sharper as the density is increased, but the one at $`r/\sigma =4`$ is not well-defined even at $`\nu =0.75`$; higher order peaks are almost invisible, indicating long-range disorder due to the shearing in contrast to the long-range order that evolves in the non-sheared system. The comparison between the sheared and non-sheared cases suggest that the structure-formation is much slower in a sheared fluid which, in turn, implies that, as expected, the freezing-point density of the former would be larger than that of the latter. The splitting of the second-peak in the non-sheared case corresponds to the onset of freezing transition34; with increasing density this splitting becomes much more prominent, with similar structural-features appearing at the successively higher-order peaks. For the sheared system, however, we do not observe similar splitting, rather we see a sharp second-peak. This also occurs for a highly inelastic sheared fluid at a much lower density33; hence, a sharp second peak is a signature of short-range ordering due to the dissipative particle-clustering. The higher-order peaks in the sheared case become prominent only if we go beyond $`\nu =\pi /4`$, which corresponds to the limit of perfect square-packing, but it is difficult to maintain homogeneous-shearing at such high densities, because the system splits into two parts, a dense, cold, crystalline area and a dilute, hot, fluid area – see below for details.

Figure 6 shows the variations of the pressure and viscosity functions ($`f_p=p/\rho T`$ and $`f_\mu =\mu /\rho \sigma \sqrt{T}`$, respectively) with density at a restitution coefficient $`e=0.7`$. We observe that both increase monotonically with density, much beyond $`\nu =0.7001`$, and we did not find the hysteretic van-der-Walls loop in our pressure data upto $`\nu =\pi /4`$, another indication that the crystallization is hindered/delayed by shear. These observations, together with our result that in the elastic limit $`𝒩_1`$ changes sign at a much lower density ($`\nu 0.62`$, see Fig. 3) than the corresponding freezing-density, suggest that the sign-reversal of $`𝒩_1`$ is not related to the freezing-transition.

#### III.2.2 Collision angle distribution

Associated with the sign-reversal of the first normal stress difference is a change in the relative magnitudes of the normal stress components ($`P_{xx}`$ and $`P_{yy}`$). Subtle changes in the direction and magnitudes of the collisional-mode of momemtum transfer could influence the individual components of the stress tensor. In order to test this hypothesis, we focus on the collision angle distribution function, $`C(\theta )`$, which is defined such that $`C(\theta )d\theta `$ is the probability of collisions occuring at an angle lying between $`\theta `$ and $`\theta +d\theta `$, with the angle $`\theta `$ being measured in the anticlockwise direction from the positive $`x`$-axis (see Fig. 13). For a fluid in equilibrium, all collisions are equally likely, and hence $`C(\theta )`$ is a uniform function of $`\theta `$, i.e. $`C(\theta )=1/\pi 0.318309`$. For a non-equilibrium system (e.g. shear flow), however, preferred collisions are dictated by the nature of the external field, leading to an anisotropic distribution for $`C(\theta )`$35-38. Following Savage & Jeffrey36 and Campbell & Brennen37, an explicit expression for $`C(\theta ;\nu ,e)`$ can be derived for the case of uniform shear flow as detailed in the Appendix. Note that the angular dependence of $`C(\theta ;\nu ,e)`$ will be modified by both the density and the restitution coefficient.

Figure 7 shows the comparison of our simulation data on $`C(\theta )`$ with the theoretical predictions of equation (29) for two different values of the restitution coefficient at a density $`\nu =0.6`$. It is observed that the probability of collisions is higher on the upstream-faces of the colliding particles, i.e. for $`\theta [\pi /2,\pi ]`$ and $`\theta [\pi /2,0]`$ (i.e. the hatched-areas on the test-particle in Fig. 13). This is a consequence of the imposed shear-field which compresses the flow-structure along the $`3\pi /4`$-direction and stretches it along the $`\pi /4`$-direction. Regarding the comparison with theory, there is, clearly, a phase-difference between theory and simulation, and the overall agreement is only qualitative.

Note in Fig. 7 that the probability of collisions on the upstream-faces increases further as the restitution coefficient increases. This, in turn, suggests that particle motion becomes more streamlined (i.e. along the streamwise direction) with increasing dissipation levels, which will naturally lead to a reduction in the transverse component of the fluctuation velocities of the particles. Thus, the macroscopic manifestation of such microscopic streamlined-motion would be an increase in the magnitude of the kinetic component of the first normal stress difference ($`𝒩_1^k`$).

Turning our attention to the range of densities where $`𝒩_1`$ undergoes a sign-reversal, we show the collision angle distributions $`C(\theta )`$ in Fig. 8 as polar plots with $`e=0.7`$; the corresponding densities are as in the subplots of Fig. 4. It is observed that the anisotropic structure of $`C(\theta )`$ gets further modified in this regime, with distinct peaks appearing near $`\theta =0`$ and $`2\pi /3`$ (see subplot $`b`$). While the peak at $`\theta =0`$ corresponds to head-on collisions between particles in the same-layer, the one at $`\theta =2\pi /3`$ clearly signals the onset of triangular-structure formation. Another noteworthy point is that the collisions on the downstream-faces of the colliding particles are rare at these densities and hence $`C(\theta )`$ can be approximated solely by its contributions from the second- and fourth-quadrants ($`\theta [\pi /2,\pi ]`$ and $`\theta [\pi /2,0]`$, respectively).

Since the momentum transfer occurs mainly due to collisions in the dense regime, the stress tensor can be approximated by

$$𝐏(𝐤𝐤)C(\theta )d\theta ,$$

where $`𝐤`$ is the unit vector joining the line of centres of the two colliding disks. Assuming now that all the collisions would occur at some average collision angle $`\theta _{av}`$ so that $`C(\theta )=C(\theta _{av})`$, and recalling that $`C(\theta )`$ is well represented in this regime by restricting $`\theta `$ only in the second- and fourth-quadrants, the expression for the first normal stress difference simplifies significantly to

$$𝒩_1\left[(𝐤𝐤)_x(𝐤𝐤)_y\right]_{\theta =\theta _{av}}C(\theta _{av}).$$
(10)
It is trivial to check that $`𝒩_1=0`$ at $`\theta _{av}=\pi /4`$. From our simulation data, we have calculated $`\theta _{av}`$ by averaging $`C(\theta )`$ over the second- and fourth-quadrants, whose variation with density is plotted in Fig. 9 for two restitution coefficients. It is observed that $`\theta _{av}`$ crosses through $`\pi /4`$ (i.e. $`3\pi /4`$) at around the critical density $`\overline{\nu }`$ for all restitution coefficients. For example, $`\theta _{av}45.16^{}`$ and $`45.04^{}`$ at $`\overline{\nu }=0.725`$ and $`0.67`$, respectively, where $`𝒩_1`$ changes sign. Thus, the microstructural signature of the sign-reversal of $`𝒩_1`$ is directly correlated with the average collision angle being greater or less than $`\pi /4`$ (or $`3\pi /4`$).

### III.3 Crystallization: Fluid-solid coexistence

Figure 10($`a`$) shows a snapshot of the system at the steady state with parameter values being set to $`\nu =0.8`$ and $`e=0.9`$. It is observed that a solid-layer coexists with two fluidized zones on either side of it. A closer look into the solid-layer reveals that the particles are arranged in a triangular-packing, representing a crystal, and thus we have a clear evidence for fluid-solid coexistence. The corresponding instantaneous streamwise velocity profile at $`t=390`$ (i.e. the image-boxes have moved $`390`$ strain units from their original position) is shown in Fig. 10($`b`$); the coarse-graining is done by binning the system into 20 equal-size bins in the transverse direction and then taking averages over all the particles in each bin. Clearly, the shearing is inhomogeneous across the Couette-gap: the crystal is aligned along the streamwise direction and hence we call it a layered-crystal; the shear-rate in the fluidized regimes on either side of the crystal is almost uniform. Note that the asymmetric nature of the velocity profile also signals the breakdown of the Lees-Edwards boundary condition as a motor for the homogeneous shear. The formation and the time-evolution of this crystal can be ascertained from Fig. 10($`c`$) which shows the corresponding evolution of the streamwise velocity at early times. We observe that the crystal has fully formed at $`t=76`$, and the velocity profile remains antisymmetric about $`y=0`$ till $`t=150`$. With further time-evolution, however, the crystal does not remain stationary in the transverse direction, rather it moves slowly with particles diffusing across the fluid-solid interface. The overall life-time of this crystal is several orders of magnitude larger than the external time-scale $`\gamma ^1`$, imposed by the shear. The corresponding collision-angle distribution $`C(\theta )`$ in Fig. 10($`d`$) shows three distinct peaks at $`\theta =0`$, $`\pi /3`$ and $`2\pi /3`$. Note that the peak at $`\theta =\pi /3`$ does not exist in the fluid phase (refer to Figs. 8$`a`$), and this provides evidence that the particles in the crystalline-phase are arranged in the triangular-packing structure. The large area of the crystalline phase as well as its relatively large value of solids fraction (hence, small mean free time) further suggest that most of the collisions occurred in the crystalline-phase.

We should remark here that, at such high densities, the inelastic collapse39-41 would eventually terminate the evolution of the system. We used the TC-model40 to avoid inelastic collapse, but could not altogether eliminate it within the crystal after some time. But the important point to note is that the crystalline-phase can be maintained for a long period of time ($`t>100`$ strain units), and hence the reported results are not transient effects.

Analogous plots for a nearly elastic system ($`e=0.99`$) are shown in Fig. 11($`a`$$`c`$) for the same density $`\nu =0.8`$. The overall features are similar to that for $`e=0.9`$, but the width of the crystalline zone is a little larger. By decreasing the dissipation-level to $`e=0.7`$, we did not observe crystal formation, with other parameters being fixed; by increasing the system-size to $`N=4096`$, however, we observed layered crystal at $`e=0.7`$. Thus, the formation of such layered-crystalline structure depends crucially on the system size and the dissipation level: the larger the system-size and/or the weaker the dissipation, the more susceptible the system is to crystallize.

Note that even if we are well below the limit of perfect square-packing ($`\nu _{sp}=\pi /60.785`$), the system could crystallize if the dissipation-levels are low; for example, we observed layered crystalline structures at $`\nu =0.75`$ with $`e=0.99`$ and $`N=4096`$. Decreasing the coefficient of restitution to $`e=0.9`$, the flow-field remained homogeneous. Thus, our layered crystalline structures appear to be tied to a long-wave instability of the elastic hard-sphere fluids56. Since our results are not driven by the inelastic dissipation, they are distinctly different from the layered shear-banding patterns, as predicted by the kinetic theory models15,42,43, in granular Couette flows. Having said that we need to mention that such dissipative layering patterns (i.e. those which become stronger with increasing dissipation levels15) were also found in simulations of a dilute sheared granular fluid43, but the simulations were allowed to evolve from an unsheared initial configuration with a special kind of boundary condition. In contrast, all our simulations started from a uniform shear condition. In the present contribution, we have mainly focussed on the non-dissipative layering in the dense limit, and the related issues of dissipative layering are relegated to a future study.

Before moving further, we make a qualitative comparison with the earlier simulation work of Campbell & Brennen37 on the bounded Couette flow of inelastic, frictional hard disks. They also reported similar layered-microstructure but due to the small system-size ($`N=40`$) and boundary-effects (they considered a shear-flow bounded by frictional walls), the distinct shear-band formation that we have reported is not evident in the snapshots of their simulations. Nevertheless, we believe that our results are akin to that reported by Campbell & Brennen. One of the referees has drawn our attention to the recent work of Campbell44 who probed the dense limit of a three-dimensional Couette shear flow using soft-sphere simulations. He was able to maintain uniform shearing even at a density of $`\nu =0.62`$, but beyond that he reported shear-band formation in that the flow-field degenerates into sheared and non-sheared zones in the gradient direction. This simply suggests that one can maintain uniform shearing in a three-dimensional geometry even beyond the analogous square-packing limit ($`\varphi =\pi /40.52`$) since the particles have now an additional degree of freedom, orthogonal to the shear plane, to rearrange themselves. Moreover, most of Campbell’s simulations were done with $`1000`$ particles (equivalently, $`100`$ particles in 2D) at a restitution coefficient of $`0.7`$. With these parameter values, we did not find layered crystals at $`\nu =0.8`$ in two-dimensions (see discussion in the next section).

#### III.3.1 Shear-induced ordering and Reynold’s dilatancy

One of the referees has drawn our attention to the recent work of Lutsko45 who studied the shear-induced ordering in a low-density ($`\nu 0.26`$) elastic hard-sphere fluid. The earliest simulations of Erpenbeck46 on 3D elastic hard-sphere fluids showed that at high shear rates the system breaks down into orderded (solid) and disordered (fluid) phases in the direction of the mean vorticity (i.e. normal to the shear-plane). This induces a long-range two-dimensional ordering, called a string-phase, in the shear-plane. More importantly, such ordering occurs only for a range of shear-rates– the lower the density, the larger this shear-rate interval. However, the latter work of Evans and Morriss47 showed that the disorder-order transition of Erpenbeck arises due to the profile-biased-thermostat47 since the string-phase vanishes completely when a profile-unbiased-thermostat47 is used. Thus, the shear-induced-ordering in elastic hard-sphere fluids appears to depend on the choice of the thermostat. Note further that all the above works probed moderately dense systems only, well below the corresponding square-packing limit. Now to compare these results with our observations in granular fluids, we first need to define an equivalent shear rate since the granular energy and the shear rate are dependent on each other for the latter system. From the energy balance equation, it is trivial to show that the granular energy has the following functional relation with the shear rate and the restitution coefficient:

$$T\frac{\gamma ^2}{1e^2}\gamma _{}^{}{}_{}{}^{2},$$

where $`\gamma ^{}=\gamma /(1e^2)^{1/2}`$ is defined as the reduced shear rate. In the quasielastic limit ($`e1`$) the reduced shear rate (and hence the granular energy) approaches infinity, which is equivalent to the high shear-rate limit of an elastic hard-sphere fluid. Since we observed layered crystalline structures only in the quasielastic limit, we may thus conclude that such structures would also persist in elastic hard-disk fluids at large shear rates.

It is interesting to ask whether the layered crystalline-structures of Figs. 10 and 11 are, in any way, related to the concept of the Reynold’s dilatancy17,18 which is explained schematically in Fig. 12. The top two sketches depict the classical constant-load48 shear-cell experiment (in two-dimensions) in which the material is subjected to a constant normal load. Fig. 12($`a`$) corresponds to an ideal situation of maximum packing, with the particles being arranged in a triangular lattice ($`\nu =\nu _{tp}=\pi /2\sqrt{3}`$). Clearly, in this situation the top and the bottom plates of the shear cell will simply slide over the material, without deforming it. However, if one of the plates is allowed to move in the vertical direction and thereby allowing the particles to rearrange themselves, e.g. as in Fig. 12($`b`$), the material can be deformed even homogeneously (if $`\nu \nu _{sp}=\pi /6`$). This is the shear-coupled volume change, commonly known as the Reynold’s dilatancy. Note, however, that our simulations mimic constant-volume experiments, since the volume of the computational box remains fixed (see Figs. 12$`c`$ and 12$`d`$); but, of course, now the pressure (i.e. normal load) can vary in response to particle motions inside the shear cell. For this case also, the uniform shearing is possible if and only if the overall density remains below the square-packing limit ($`\nu <\nu _{sp}`$). However, for densities above this value ($`\nu _{sp}<\nu <\nu _{tp}`$), the shearing can be started only if we allow the particles to rearrange themselves. This is possible if a part of the system becomes denser, allowing free volumes to the rest of the system which is nothing but Reynold’s dilatancy too. Hence, we will end up with a crystalline-phase coexisting with a fluid phase, a typical example of which is shown schematically in Fig. 12($`d`$). Thus the phenomenon of Reynold’s dilatancy, for densities $`\nu _{sp}<\nu <\nu _{tp}`$, would make the ordering transition, as depicted in Figs. 10 and 11, more prominent, and this effect would be much stronger in two-dimensions.

## IV Consequences for the constitutive modeling: Relaxation models

Here we attempt to describe the normal-stress behaviour of a granular fluid using the standard relaxation-type models. Prior literature on the dense-gas kinetic theory, which forms the foundation of theoretical developments of granular fluids in the rapid-shear regime, indicates that such a stress relaxation mechanism does also exist in granular fluids49-51. The relaxation-type models are routinely used to describe the non-Newtonian behaviour of viscoelastic/viscoplastic materials, and hence might be apt for granular fluids in the dense limit as well. The recent work of Zhang & Rauenzahn52-53 clearly shows that such viscoelastic stress relaxation mechanism exists in dense granular flows. Following a rigorous statistical mechanical procedure, they derived an evolution equation for the collisional stress tensor which boils down to a frame-indifferent viscoelastic model, with the Jaumann derivative appearing directly without appealing to objectivity arguments.

Let us consider the viscoelastic relaxation approximation suggested recently by Jin & Slemrod51 to regularize the Burnett order equations of Sela & Goldhirsch5 for a low-density granular fluid. Their proposed equation for the pressure deviator, in our notation, is

$`𝚷+\tau _1\left({\displaystyle \frac{D𝚷}{Dt}}𝐋^T𝚷𝚷𝐋+\frac{2}{d}tr(𝚷𝐋)\mathrm{𝟏}\right)`$
$`+\tau _2\left(𝐒𝚷+𝚷𝐒\frac{2}{d}tr(𝚷𝐒)\mathrm{𝟏}\right)=𝚷^{eq}`$ (11)
where

$`𝚷^{eq}`$ $`=`$ $`2\mu 𝐒\left(\lambda 𝐮\right)\mathrm{𝟏}+𝚷_2+𝚷_3,`$
$`𝐒`$ $`=`$ $`\frac{1}{2}\left(𝐋+𝐋^T\right)\frac{1}{d}\left(𝐮\right)\mathrm{𝟏},`$
$`𝐋`$ $`=`$ $`(𝐮)^T=\left({\displaystyle \frac{u_i}{x_j}}\right),`$
$`\lambda `$ $`=`$ $`(\zeta {\displaystyle \frac{2}{d}}\mu ),`$
$`\tau _1`$ $`=`$ $`0.3211\left({\displaystyle \frac{\mu }{p}}\right),`$
$`\tau _2`$ $`=`$ $`0.58775\left({\displaystyle \frac{\mu }{p}}\right).`$
Here $`\tau _1`$ and $`\tau _2`$ are relaxation times, $`d`$ the dimensionality of the system, $`𝐋`$ is the velocity gradient, $`𝐒`$ the deviatoric part of the rate of strain tensor, $`\mu `$ the shear viscosity, $`\zeta `$ the bulk viscosity and $`\mathrm{𝟏}`$ the identity tensor; $`𝚷_2`$ and $`𝚷_3`$ are higher order terms, explicitly written down in Jin & Slemrod51. Note that both relaxation times are proportional to the ratio of the shear viscosity and the pressure, and hence proportional to the mean free time. In the limits of $`\tau _1`$, $`\tau _20`$ and $`𝚷_2`$, $`𝚷_30`$, we recover the standard Newtonian model for the stress tensor.

Neglecting the higher-order terms, an expression for the first normal stress difference can be obtained for the steady uniform shear flow:

$$𝒩_1=\left[\frac{4\tau _1}{2(1+\tau _1^2)+\tau _2(1+\tau _1\tau _2)}\right]\left(\frac{\mu }{p}\right)>\mathrm{\hspace{0.33em}\hspace{0.33em}0}.$$
(12)
This quantity is always positive, as in our simulation results for dilute flows, if the two relaxation times are of the same order.

It is important to note that the above evolution equation does not satisfy the principle of material frame indifference (MFI) which states that the constitutive laws should be invariant under rigid-rotation7,8. The scalar field $`\varphi `$, the vector field $`𝐯`$ and the tensor field $`𝚷`$ are called frame-indifferent or objective if the following relations hold for all $`t`$:

$`\varphi ^{}(𝐱^{},t^{})`$ $`=`$ $`\varphi (𝐱,t),`$ (13)
$`𝐯^{}(𝐱^{},t^{})`$ $`=`$ $`𝐐(t)𝐯(𝐱,t),`$ (14)
$`𝚷^{}(𝐱^{},t^{})`$ $`=`$ $`𝐐(t)𝚷(𝐱,t)𝐐(t)^T,`$ (15)
where $`\varphi `$, $`𝐯`$ & $`𝚷`$ and $`\varphi ^{}`$, $`𝐯^{}`$ & $`𝚷^{}`$ are defined in two different frames $``$ and $`^{}`$, respectively, and $`𝐐(t)`$ is a proper orthogonal tensor. Here $`^{}()`$, with $`()`$ denoting the set of all frames obtainable from a given frame $``$ by observer transformations. That the stress-tensor in a granular gas is not a frame-independent quantity (as in the hard-sphere gas54,55) is well-known. Since the kinetic component of the first normal stress difference remains positive at all densities, the kinetic stress tensor can be modelled using a non-objective equation as discussed above. For the collisional stress tensor, one can postulate a similar evolution equation as in Eqn. (11), but $`\tau _1`$ must be multiplied by a factor which must change sign at the critical density. However, this equation would remain frame-dependent even in the dense limit. Thus, if one has to recover the Boltzmann-limit of relaxation-type equations, a frame-indifferent approach does not appear to work.

A question now arises as to the possibility of modelling normal stress differences using the standard frame-indifferent relaxation type models. The simplest way is to use either the lower-convected or the upper-convected equations for the pressure deviator:

$`𝚷+\alpha (\nu ,e)\tau \left({\displaystyle \frac{\mathrm{D}𝚷}{\mathrm{D}t}}+\left\{𝐋^T𝚷+𝚷𝐋\frac{2}{d}\mathrm{tr}(𝚷𝐋)\mathrm{𝟏}\right\}\right)`$ $`=`$ $`2\mu 𝐒\left(\lambda 𝐮\right)\mathrm{𝟏},`$ (16)
$`𝚷+\alpha (\nu ,e)\tau \left({\displaystyle \frac{\mathrm{D}𝚷}{\mathrm{D}t}}\left\{𝐋𝚷+𝚷𝐋^T\frac{2}{d}\mathrm{tr}(𝚷𝐋)\mathrm{𝟏}\right\}\right)`$ $`=`$ $`2\mu 𝐒\left(\lambda 𝐮\right)\mathrm{𝟏},`$ (17)
respectively. Here $`\alpha (\nu ,e)`$ is an empirical constant, dependent on both the density and restitution coefficient. For both cases, the first normal stress difference is

$$𝒩_1=\frac{2\tau \alpha }{1+\tau ^2\alpha ^2}\left(\frac{\mu }{p}\right).$$
(18)
Clearly, if $`\alpha (\nu ,e)`$ is obtained from simulation, its sign-reversal would also correspond to the sign-reversal of $`𝒩_1`$. Comparing the Jin-Slemrod equation with its corresponding frame-indifferent analog (lower-convected model), we conclude that the loss of frame-indifference shows up as a sign-change of the first normal stress difference. It would be interesting to investigate whether one could relax the Sela-Goldhirsch equations4 using a frame-indifferent approach without violating the entropy inequality51.

Similarly, one could postulate evolution equations using other objective derivatives. In this regard, the co-rotational Jeffrey’s model seems to be the ideal choice:

$$𝚷+\tau _1(\nu ,e)\frac{𝒟𝚷}{𝒟t}=2\mu \left(𝐒+\frac{\lambda }{2\mu }\left(𝐮\right)\mathrm{𝟏}+\tau _2(\nu ,e)\frac{𝒟𝐒}{𝒟t}\right)$$
(19)
with $`𝒟/𝒟t`$ being the Jaumann derivative7. The corresponding first normal stress difference is

$$𝒩_1=\frac{2(\tau _1\tau _2)}{1+\tau _1^2}\left(\frac{\mu }{p}\right),$$
(20)
which is positive/negative depending on whether $`\tau _1`$ is less/greater than $`\tau _2`$. Thus, the frame-indifferent relaxation models are able to predict positive and negative first normal stress differences. (For the steady homogeneous shear flow, one can also model positive/negative normal stress differences by postulating a general orthonormal basis, generated by the nilpotent basis tensors, which satisfies the objectivity requirement; for related issues, the reader is referred to Goddard8.)

## V Summary and Conclusion

We have probed the non-Newtonian behaviour and the incipient crystalline-phase in a sheared, monodisperse, two-dimensional granular fluid. The standard event-driven technique is used to simulate a box of hard-disks under homogeneous shear deformation. The information about the stress tensor is obtained by decomposing it in the standard way: $`𝐏=p\mathrm{𝟏}+𝚷`$, where $`p`$ is the pressure and $`𝚷`$ the pressure deviator. The non-Newtonian behaviour is quantified in terms of the first normal stress difference: $`𝒩_1=(\mathrm{\Pi }_{xx}\mathrm{\Pi }_{yy})/p`$.

The granular fluid is non-Newtonian with a measurable first normal stress difference ($`𝒩_1`$) which is positive (if the stress is defined in the compressive sense) in the dilute limit. Interestingly, however, $`𝒩_1`$ changes from positive to negative at a critical density in the dense regime. By decomposing $`𝒩_1`$ into the kinetic and collisional contributions, $`𝒩_1=𝒩_1^k+𝒩_1^c`$, we found that while $`𝒩_1^k`$ is always positive and decays to zero in the dense limit, $`𝒩_1^c`$ has a non-monotonic variation with density. In particular, $`𝒩_1^c`$ increases from zero in the dilute limit as $`\nu `$ increases, reaches a maximum at some value of $`\nu `$ and then decreases, eventually becoming negative in the dense limit. The density at which $`𝒩_1^c=0`$ ($`\nu \overline{\nu }`$) depends crucially on the level of micro-scale dissipation; in particular, $`\overline{\nu }`$ increases as the coefficient of restitution decreases. We have constructed a phase-diagram in the ($`\nu ,e`$)-plane by identifying the regions where $`𝒩_1`$ is positive/negative.

We have shown that the origin of the first normal stress difference, in the dense limit, is tied to shear-induced collisional anisotropies. The underlying mechanism is distinctly different from that is known for a dilute granular gas4,16 where the anisotropy in the second moment of the fluctuation velcoity, which is a Burnett-order effect, gives rise to normal stress differences.

At the micro-level, the particles undergo reorganization as the dense-limit is approached. The signatures of microstructural-reorganization have been captured by probing the collision-angle distribution, $`C(\theta )`$, which is anisotropic due to the presence of the mean shear. In particular, we have found that the topology of the anisotropic-structure of $`C(\theta )`$ changes, with collisions occurring at certain preferred angles on the upstream-faces of the colliding pairs. The sign-reversal of $`𝒩_1`$ is correlated with a preferred value of the average collision angle, $`\theta _{av}=\pi /4\pm \pi /2`$, averaged over the upstream-faces of the colliding particles.

The time evolution of the sheared granular fluid leads to crystallization in the dense limit, signalling the coexistence of fluid and solid. The particles are arranged in a triangular-packing inside the crystal, and it moves as a layer in the streamwise direction. The formation of such layered-crystalline structure depends crucially on the system size and the dissipation level: the larger the system-size and the weaker the dissipation, the more susceptible the system is to crystallize. This appears to be related to a long-wave instability56 of the elastic hard-sphere fluids.

The present work clearly shows that the available kinetic-theory-based rheological models for granular fluids are not adequate to predict the behaviour of the first normal stress difference in the dense limit. Certain microstructural-features, like the preferred distribution of collisions which eventually leads to crystal-formation, should be incorporated into the theory. At such high densities, many-body effects (both positional and velocity correlations) are important57,41 and the appropriate kinetic description is the BBGKY-hierarchy57. To incorporate the observed normal stress behaviour into the framework of plausible constitutive models, we showed that the standard frame-indifferent relaxation type models can be used to model both positive and negative first normal stress differences. In this regard, the two-parameter Jeffrey’s model7 appears to be the ideal choice; however, we are unable to recover the corresponding Boltzmann limit which is known to be non-objective. On the whole, we believe that a lot remains to be done for a better understanding of the dense-phase rheology of granular fluids even in the hard-sphere limit.

## VI Acknowledgement

M.A. acknowledges the financial support by the Alexander von Humboldt Foundation and S.L. acknowledges the support of the Deutsche Forschungsgemeinschaft. M.A. is grateful to Hans J. Herrmann for discussions and hospitality at ICA1, and also acknowledges discussions with Joe D. Goddard on certain aspects of normal stress behaviour in granular fluids. We thank John F. Brady for directing us to relevant references on normal stresses in suspensions. We also thank three anonymous referees for their constructive comments.

## Appendix A Shear-induced anisotropy and the collision angle distribution

Here we are interested in the shear-induced anisotropy of the collision angle distribution of an inelastic hard-disk fluid. Following Savage & Jeffrey36, an expression for the collision angle distribution $`C(\theta )`$ is derived, which is compared with the simulation data in IIB.2.

To calculate the probablility of collisions at a specific angle $`\theta `$, we focus on Fig. 13 with two particles colliding at $`𝐫`$. Note that $`\theta `$ is measured anticlockwise from the positive $`x`$-axis. For collisions to occur in a time $`\delta t`$, the center of particle $`1`$ must lie inside the volume $`\sigma \delta 𝐤(𝐪𝐤)\delta t`$, where $`𝐪=𝐜_1𝐜_2`$ is the relative velocity of the colliding pair. Thus the expected total number of collisions (per unit time and unit area) with the line of centres $`𝐤`$ lying between $`𝐤\delta 𝐤/2`$ and $`𝐤+\delta 𝐤/2`$ is given by

$$\sigma f^{(2)}(𝐜_1,𝐫_1,𝐜_2,𝐫_2)(𝐪𝐤)(𝐤𝐧)𝑑𝐤𝑑𝐜_1𝑑𝐜_2,$$
(21)
where $`f^{(2)}()`$ is the two-particle distribution function which is defined so that $`f^{(2)}(𝐜_1,𝐫_1,𝐜_2,𝐫_2)`$ $`d𝐜_1d𝐜_2d𝐫_1d𝐫_2`$ is the number of pairs of particles such that the particle $`i`$ is located in an area element $`d𝐫_1`$ about $`𝐫_1`$ with its velocity in the interval $`d𝐜_1`$ about $`𝐜_1`$ while particle $`j`$ is located in an area element $`d𝐫_2`$ about $`𝐫_2`$ with its velocity in the interval $`d𝐜_2`$ about $`𝐜_2`$. To progress further, we have to invoke the assumption of molecular chaos and hence the expected number of collisions is

$$\sigma g(𝐫_1,𝐫_2)f^{(1)}(𝐜_1,𝐫_1;𝐮(𝐫_1))f^{(1)}(𝐜_2,𝐫_2;𝐮(𝐫_2))(𝐪𝐤)(𝐤𝐧)𝑑𝐤𝑑𝐜_1𝑑𝐜_2,$$
(22)
where $`g(𝐫_1,𝐫_2)`$ is the pair-distribution function. For the steady uniform shear flow, $`g(𝐫_1,𝐫_2)`$ is calculated from the relation35,36,58:

$$g(𝐫_1,𝐫_2)=\frac{2g_c(\nu )}{n^2}_{𝐪𝐤>0}f^{(1)}(𝐜_1,𝐫_1;𝐮(𝐫_1))f^{(1)}(𝐜_2,𝐫_2;𝐮(𝐫_2))𝑑𝐜_1𝑑𝐜_2,$$
(23)
where $`g_c(\nu )`$ is the contact value of the pair-distribution function and $`𝐪𝐤>0`$ implies that the integration be carried out for impending collisions.

As a first approximation, the single particle velocity distribution function $`f^{(1)}(𝐜_1,𝐫_1;𝐮(𝐫_1))`$ is assumed to have the Maxwellian-form:

$$f^{(1)}(𝐜_1,𝐫_1;𝐮(𝐫_1))=\left(\frac{nm}{\pi k_BT}\right)\mathrm{exp}\left[\frac{m\left(𝐜_1𝐮(𝐫_1)\right)^2}{2k_BT}\right],$$
(24)
where $`T`$ is the granular temperature (fluctuating kinetic energy) and $`k_B`$ the usual Boltzmann constant. Now transforming the particle velocities ($`𝐜_1`$, $`𝐜_2`$) to their center-of-mass and relative velocities, equation (23) can be integrated to arrive at the following explicit expression for the pair-distribution function36:

$$g(𝐫_1,𝐫_2)=g_c(\nu )\mathrm{erfc}\left[\frac{2𝐤𝐮(𝐫_2)}{(2k_BT)^{1/2}}\right],$$
(25)
where $`\mathrm{erfc}()`$ is the complementary error function. Using the above expression for the pair-distribution function and transforming in terms of polar coordinates ($`r,\theta `$), the integral for the normalized collision angle distribution yields37

$$C(\theta )=𝒜(T)\left[\mathrm{exp}\left(\frac{\mathrm{sin}^2\theta \mathrm{cos}^2\theta }{2T}\right)\frac{g(\theta )\mathrm{sin}\theta \mathrm{cos}\theta }{\sqrt{T}}\right]g(\theta ),$$
(26)
where $`g(\theta )`$ is the angular pair-distribution function given by

$$g(\theta )\frac{g(𝐫_1,𝐫_2)}{g_c(\nu )}=\mathrm{erfc}\left[\frac{\mathrm{sin}\theta \mathrm{cos}\theta }{\sqrt{2T}}\right]$$
(27)
and $`𝒜(T)`$ is a normalization constant.

For the uniform shear flow, an expression for the granular temperature, can be obtained from the energy balance equation, by equating the energy production due to shear-work with the energy loss due to collisions:

$`\mu \left({\displaystyle \frac{du}{dy}}\right)^2`$ $`=`$ $`𝒟`$
$`T`$ $`=`$ $`f_\mu (\nu )/f_𝒟(\nu ,e)`$ (28)
where $`\mu =\rho _p\sigma f_\mu (\nu )\sqrt{T}`$ is the shear viscosity and $`𝒟=(\rho _p/\sigma )f_𝒟(\nu ,e)T^{3/2}`$ the collisional dissipation rate, with

$`f_\mu (\nu )`$ $`=`$ $`{\displaystyle \frac{\sqrt{\pi }\nu }{8}}\left[{\displaystyle \frac{1}{\nu g_c(\nu )}}+2+\nu g_c(\nu )\left(1+{\displaystyle \frac{8}{\pi }}\right)\right],`$
$`f_𝒟(\nu ,e)`$ $`=`$ $`{\displaystyle \frac{4}{\sqrt{\pi }}}(1e^2)\nu ^2g_c(\nu ).`$
Substituting this expression for $`T`$, the normalized collision distribution function becomes

$$C(\theta ;\nu ,e)=𝒜(T)\left[\mathrm{exp}\left(\frac{f_𝒟(\nu ,e)\mathrm{sin}^2\theta \mathrm{cos}^2\theta }{2f_\mu (\nu )}\right)g(\theta )\mathrm{sin}\theta \mathrm{cos}\theta \sqrt{\frac{f_𝒟(\nu ,e)}{f_\mu (\nu )}}\right]g(\theta )$$
(29)
and the angular pair-distribution function

$$g(\theta ;\nu ,e)=\mathrm{erfc}\left[\mathrm{sin}\theta \mathrm{cos}\theta \sqrt{\frac{f_𝒟(\nu ,e)}{2f_\mu (\nu )}}\right].$$
(30)
It is clear that the angular dependence of $`C(\theta ;\nu ,e)`$ is modified by both the inelastic dissipation and the density.