# Self-assembly and non-equilibrium phase coexistence in a binary granular mixture

Andrea Plati, Raphaël Maire, François Boulogne, Frédéric Restagno, Frank

Smallenburg, Giuseppe Foffi

► To cite this version:

Andrea Plati, Raphaël Maire, François Boulogne, Frédéric Restagno, Frank Smallenburg, et al. Self-assembly and non-equilibrium phase coexistence in a binary granular mixture. The Journal of Chemical Physics, 2025, 163 (5), pp.054509. 10.1063/5.0268711. hal-04859858

HAL Id: hal-04859858

https://hal.science/hal-04859858v1

Submitted on 5 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Self-assembly and non-equilibrium phase coexistence in a binary granular mixture

A. Plati,¹ R. Maire,¹ F. Boulogne,¹ F. Restagno,¹ F. Smallenburg,¹ and G. Foffi¹

¹Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France

(Dated: October 30, 2024)

We report the experimental observation of a square crystalline phase in a vibrated binary mixture of spherical grains. This structure spontaneously forms from a disordered state, consistently with predictions obtained in an equilibrium system with similar geometrical properties under conservative dynamics. By varying the area fraction, we also observe stable coexistence between a granular fluid and an isolated square crystal. Using realistic simulations based on the discrete element method and an idealized collisional model integrated via event-driven molecular dynamics, we not only reproduce experimental results but also help to gain further insights into the non-equilibrium phase coexistence. Through the direct phase coexistence method, we demonstrate that the system shows behavior highly similar to an equilibrium first-order phase transition. However, the crystal remains at a higher granular temperature than the fluid, which is a striking non-equilibrium effect. Through qualitative arguments and supported by kinetic theory, we elucidate the role of the coupling between local structure and energy transfer mechanisms in sustaining kinetic temperature gradients across the fluid-solid interface.

## I. INTRODUCTION

The study of spontaneous self-assembly of ordered structures starting from disordered phases unites the interests of a diverse scientific community ranging from physicists and chemists to engineers [1–6]. The paradigmatic example of hard-sphere crystallization demonstrates how equilibrium statistical mechanics can provide a general theoretical framework to understand self-assembly processes in realistic systems such as colloidal suspensions [7]. Following equilibrium thermodynamics, the spontaneous tendency of an amorphous system to form ordered structures can be understood in terms of phase stability: In a certain region of the control parameter space, ordered macrostates are more thermodynamically stable than disordered ones (i.e. they minimize the system’s free energy). In many cases, this results in a first-order phase transition which also implies a coexistence region between the two phases. Remarkably, self-assembly of ordered structures and phase coexistence have been also observed in systems that are far from thermodynamic equilibrium such as active [8–14] and granular matter [15–21]. These are intrinsically open systems where energy constantly flows in and out. Thus, despite the possibility for them to attain the so-called non-equilibrium stationary states, where the statistical properties of macroscopic observables do not change with time, any notion of thermodynamic stability based on free energy arguments is lost. Nonetheless, out-of-equilibrium self-assembly processes can be extremely similar to equilibrium ones in geometrically equivalent conditions. A striking example of that is provided by the crystallization of quasi-2D vibrated granular matter [19, 22–25]. It has been known for more than two decades that a monolayer of monodisperse spherical grains vibrated on a substrate can self-assemble into a hexagonal crystal [19, 22], and recently it was observed that bidisperse mixtures can spontaneously form a quasi-crystalline structure [25]. In these examples, the condi-

tions for the liquid-solid transition and the evolution of the structure resemble those of equivalent 2D hard-disk systems undergoing conservative dynamics. Nevertheless, dynamical observables show clear deviations from equilibrium phenomenology, such as non-Maxwellian velocity distributions [25, 26]. Additionally, the behavior of the granular crystallization turns out to be sensitive to variations in the frequency and the amplitude of the vibration, or the material properties of the grains. Relevant examples are: i) the observation of a non-equilibrium phase coexistence between a “hot” granular fluid and a “frozen” crystal [15, 27], ii) a change from a two-step liquid-solid phase transition involving an intermediate hexatic phase to a one-step first-order-like scenario [28] and iii) the suppression of ordering induced by strong inelasticity [17].

Based on these observations, it is clear that we cannot simply conclude that granular quasi-2D systems behave in the same way as equilibrium hard-disks, nor can we entirely dismiss the possibility of predicting their behavior from equilibrium arguments. Given the importance of granular matter for many applications [29–35] and the lack of a general theoretical framework for its description, understanding which physical conditions are suitable to be treated with effective equilibrium-like theories represents an intriguing open problem. In this paper, building on our earlier work on granular quasicrystals [25], we tackle this problem by studying experimentally and numerically a vibrated granular system able to self-assemble into a periodic binary crystal with square symmetry. Specifically, we demonstrate self-assembly of square crystals and stable fluid-crystal coexistences in both experiments and simulations of binary mixtures of granular beads. In the realm of granular materials, one usually expects polydisperse mixtures to undergo size segregation effects [36–41] while any crystalline structure built upon more than one particle species imply a homogeneous concentration of grains with different sizes. In light of this, a promising by-product of our study is the
application of self-assembly of binary crystals as a way
to realize homogeneous granular mixtures [42] even in
the presence of strong non-equilibrium effects that would
promote size segregation. To our knowledge, the only
previous example of periodic binary granular crystals
obtained in experiments relies on the presence of elec-
trostatic interactions between the grains [43], while the
dynamics of our system are purely collision-driven. We
will also see that our system is particularly well suited to
the study of non-equilibrium phase coexistence, a sub-
ject which has attracted considerable interest in recent
years thanks to a series of experimental and numerical
observations [9, 10, 15, 16, 20]. A striking feature of non-
equilibrium phase coexistence is the possibility of main-
taining stable kinetic temperature gradients between the
two different phases. This suggests the presence of a net
heat-flux which is not compatible with thermodynamic
equilibrium. Kinetic temperature differences in gas-solid
and liquid-solid coexistence have been observed in gran-
ular [15, 16] and active [9, 10] matter. Here, equilibrium-
like general principles able to predict their occurrence
and describe their characteristics are absent: Each spe-
cific case needs to be analyzed to understand the key
physical ingredients that make the coexistence possible
and determine its properties such as, for example, which
of the two phases is hotter. Through a numerical anal-
ysis, we show that in our systems the crystalline phase
has a higher granular temperature than the coexisting
liquid one. We explain this observation by using the rela-
tionship between the system’s local structure and energy
transfer mechanisms. Remarkably, we find that the ki-
netic theory of binary granular mixtures is able to predict
this phenomenon in a simplified granular model.

The remainder of this paper is structured as follows.
In Sec. II, we discuss our experimental results about the
self-assembly of the granular binary crystal and its stable
coexistence with a granular liquid phase. In Sec. III, we
focus on numerical results obtained following two com-
plementary approaches: realistic simulations based on
the Discrete Element Methods (DEM) and Event-Driven
Molecular Dynamics (EDMD) of a simplified 2D coarse-
grained granular model. Using these simulations, we first
confirm and extend the study on self-assembly, and then,
also with the help of kinetic theory, we perform a system-
atic study of the granular temperature in the coexistence
region. In Sec. IV we draw our conclusions and outline
possible future developments of the present study.

II. EXPERIMENTS

A. Experimental setup and preliminary study

Our experimental system (see Fig. 1) consists of a bi-
nary mixture of spherical polyamide beads (PROLABO)
composed of N<sub>S</sub> grains with diameter σ<sub>S</sub> = 2 mm and
N<sub>L</sub> with diameter σ<sub>L</sub> = 4 mm. Polyamide has a density
ρ = 1.14 g/cm³. The system is confined in a quasi-2D

square box (height *h* and width *L* ≫ *h*) which is ver-
tically vibrated by an electrodynamic shaker (Brüel &
Kjær, LDS V400). The bottom horizontal surface of the
box (1 cm thick) is made of anodized aluminium and,
as suggested by a previous study [26], has been sand-
blasted (glass sand, granulometry [425-600] µm) before
anodization to improve grain mobility in the horizontal
plane. The upper plate (2 mm thick) is made of glass
which allows for the direct visualization of the 2D pro-
jection of the system. Snapshots of the horizontal spatial
configurations can be acquired with a high-resolution (5
MP) camera (Basler a2A2590) placed at the top of the
setup. Unless differently stated, we take a snapshot every
30 seconds. The camera is equipped with a lens (Basler
C125-1218-5M-P) having a fixed focal length of 12.0 mm.
A signal generator (Keysight 33220A) is used to design a
desired waveform *a*<sub>*z*</sub>(*t*) which fixes the vertical accelera-
tion of the box (more details on the experimental setup
and procedure are given in the caption of Fig. 1, in the
Supplementary Material (SM) of this paper and in the
Supplementary Information of Ref. 25).

To characterize the geometrical properties of the mix-
ture, we define the size ratio *q* = *σ*<sub>S</sub>/*σ*<sub>L</sub>, the fraction of
small grains *x*<sub>S</sub> = *N*<sub>S</sub>/*(N*<sub>S</sub> + *N*<sub>L</sub>)* and the area fraction
*φ* = (*N*<sub>S</sub>*σ*<sub>S</sub><sup>2</sup> + *N*<sub>L</sub>*σ*<sub>L</sub><sup>2</sup>)*π/4*L<sup>2</sup>*. Once a granular mixture

<figure><img src="image_2.png" /></figure>

FIG. 1. Sketch of the experimental setup. A binary mixture of polyamide grains lies on a sandblasted anodized aluminium plate, which is vertically vibrated. Four aluminium spacers are screwed on the plate to confine the motion of the grains in a square area of side $L = 15$ cm through lateral walls of height $h = 6$ mm while the full lateral side of the plate is $L_p = 20$ cm. The granular system is confined from above by a glass cover which allows for the direct visualization of the system through a camera. To ease the imaging of the grains, four light panels are placed around the apparatus. The measurements of the plate acceleration are made possible by a one-axis accelerometer (Brüel & Kjær, Type 4534-B-001) that can be rigidly joined to its bottom side. The base of the apparatus is held up by three leveling feet (not shown), which allow for horizontal adjustment.
with given geometrical properties {$q, x_S, \phi$} is prepared, the experimental system is initialized by pouring the mixture of grains into the box and randomizing the initial configuration by spreading them manually. Then, the box is vertically accelerated following a prescribed sinusoidal signal $a_z(t) = g\Gamma\sin(2\pi ft)$. Vibrations are then characterized by the adimensional maximum acceleration $\Gamma$ and the frequency $f$ ($g$ is the gravity acceleration). In the SM, we provide more details about the tuning of the driving signals.

Our granular mixtures have been prepared in a range of $q, x_S$ and $\phi$ for which a geometrically equivalent system undergoing thermal motion exhibits the spontaneous formation of a binary crystalline phase with a squared symmetry [44, 45]. In Fig 2, we provide an essential sketch of this equilibrium reference system which consists of a binary mixture of non-additive hard disks effectively describing a system of spheres laying on a substrate. We also show a portion of the crystal self-assembled under thermal motion. For the rest of the paper, we will refer to this type of binary square crystal as S1 [44–46]. The adimensional values {$q, x_S, \phi$}, are referred to as the equilibrium or geometrical parameters. In contrast, the physical properties of the system that are not present in the hard-disk reference system (i.e. grain material, driving signals, roughness of the bottom plate, etc.) are designated as non-equilibrium parameters. In order to narrow down our parameter space, we first calibrated our driving conditions ($\{\Gamma, f\}$) by performing short self-assembly experiments with different drivings and optimizing for conditions where horizontal movement of the grains is fast, but small particles avoid stacking on top of each other. The best results were obtained with sinusoidal shaking at $\Gamma = 1.79$ and $f = 53$ Hz. Hence, we kept these values fixed for all other experimental runs.

## B. Self-assembly of the square binary crystal

The first goal of the experiment was to investigate whether it is possible to obtain the S1 phase at the macroscopic scale using non-equilibrium athermal driving (i.e. mechanical vibration). In the equilibrium reference system, S1 self-assembly has been observed for $q \in [0.45, 0.5]$, $x_S \in [0.45, 0.6]$ and $\phi \ge 0.84$ [45]. To remain in the same range of geometrical parameters, we prepared a mixture of polyamide beads with $q = 0.5$, $x_S = 0.51$, $\phi = 0.858$ and then performed a very long experiment (~90 hours) with the previously found optimal driving sinusoidal driving ($\Gamma = 1.79, f = 53$ Hz). Since size segregation phenomena are ubiquitous in granular matter [36–38, 40, 41, 47–50], we expect the self-assembly of S1 crystals (where the two species are completely mixed) to compete with the formation of hexagonal domains of small and large particles. In view of this, we characterize the evolution of our system with the square ($q_4$) and hexagonal ($q_6$) bond orientational order parameters (BOOPs). These are scalar observables ac-

<figure><img src="image_3.png" /><figcaption>FIG. 2. Sketch of the equilibrium reference system of non-additive hard disks in which the S1 phase was first observed [44, 45]. On the left-hand side, we show the mapping between spheres laying on a plane and non-additive disks which undergo elastic collision as they reach a distance $\sigma_{LS} = \sqrt{\sigma_L \sigma_S}$. On the right-hand side, we show a portion of the S1 crystal assembled in the equilibrium reference system.</figcaption></figure>

counting for square and hexagonal crystallinity defined as:

$$q_k = \frac{1}{N_L} \sum_{m=1}^{N_L} q_k^m, \quad q_k^m = \frac{1}{N_m} \left| \sum_{j=1}^{N_m} e^{ik\theta_{mj}} \right| \quad (1)$$

with $k = \{4, 6\}$ and where $N_m$ is the number of nearest neighbours (identified through a cutoff distance) of the $m$th grain and $\theta_{mj}$ denotes the angle of the vector joining the centers of two particles with respect to a reference direction. In Eq. (1), we consider only large particles as these are sufficient for the characterization of the local orientational order. To focus specifically on the S1 crystalline phase, we also introduce a different scalar observable namely the fraction $F_{S1}$ of large particles belonging to an S1 environment. This is obtained by counting the particles with $q_4^m > q_4^*$ having at least three nearest neighbours with $q_4^m > q_4^*$, where $q_4^*$ is an arbitrary threshold that we chose equal to 0.8. $F_{S1}$ is more precise than the global $q_4$ as an indicator of crystallization since it does not take into account the contribution of isolated particles with a high $q_4^m$ (which can also be present in a fluid phase). In Fig. 3a, we plot $F_{S1}$ as a function of time and show a zoom of an S1 unit cell formed in the experiment; in Fig. 3b-e, we provide the direct visualization of the self-assembly process at four consecutive times. At the beginning, the system is completely amorphous with $F_{S1} \sim 0$ (Fig. 3b), then, for the first 40 hours, we observe a time evolution of $F_{S1}$ consisting of long stationary periods interrupted by relatively fast increases which correspond to the growth of large S1 clusters with different orientations (Fig. 3c and 3d). At the same time, the BOOPs plotted in the bottom right inset show an increase in the
<figure><img src="image_4.png" /><figcaption>FIG. 3. Results from a long experiment with polyamide beads on a rough surface under sinusoidal vibration. Here $q = 0.5$, $x_S = 0.51$, $\phi = 0.858$, $\Gamma = 1.79$, $f = 53$ Hz. a) Time evolution of the fraction of particles belonging to an S1 environment (see the text for the definition). The red dots mark the time at which the snapshots of the system shown in panels b-e) are taken (they are in chronological order). In the upper-left inset of panel a), we show the time evolution of the average large particle displacement over a time interval of 30 seconds (smoothed using a running average with a time window of 5 minutes). We use this quantity to characterize the average mobility of the system. In the bottom-right inset of panel a), we plot the time evolution of $q_4$ and $q_6$. Also shown here is a zoom-in on the unit cell of the S1 phase, which consists of four large grains surrounding a small one.</figcaption></figure>

average local square orientational order and a decrease in the hexagonal one. We point out that the main S1 forma- tion is observed around the few largest clusters nucleated at relatively short times. The remaining time of the ex- periment is characterized by a very slow increase of $F_{S1}$ during which the two largest clusters coarsen into a single S1 domain spanning a macroscopic fraction of the system (Fig. 3e). We also observe the formation and growth of smaller S1 clusters in the bottom part of the system. As expected, the self-assembly process is accompanied by a gradual decrease in particle mobility, which we define as the average absolute displacement of large particles over a time interval of 30 seconds (upper-left inset of Fig. 3a). During the second half of the experiment, we note the ap- pearance of hexagonal domains of small and large grains in different parts of the system. This is reflected also in the time evolution of the BOOPs where we see an in- crease of $q_6$ starting around $t = 60$ h. We can explain this effect by considering that the S1 structure requires a homogeneous composition of small and large grains; therefore, starting from a local amorphous environment with a non-perfectly homogeneous composition, the for- mation of S1 will tend to expel the excess species. This in

turn can lead to local regions with a high concentration of one of the species, which favours the formation of hexag- onal domains. On top of that, over very long experiments the horizontal calibration of the system may slowly get lost, introducing the effect of gravity and inhomogeneity in the granular dynamics. This promotes size segregation [37, 38, 41] and may therefore play a role in the formation of hexagonal domains. However, the important thing to note here is that these effects are rather marginal com- pared to the mixing of the two species within the S1 structure which continues to grow despite the formation of competing structures. The reproducibility of this re- sult is supported by a shorter experiment (30 hours) with a slightly lower area fraction ($\phi = 0.851$) that exhibited the formation of S1 domains in the early dynamics (see SM).

Finally, we remark that mixing particles of different sizes in a granular system is a paradigmatic problem par- ticularly relevant in different applications [37, 39, 51]. As argued in a recent numerical study [42], the self-assembly of binary crystals in a polydispersed granular system can be a way to address the problem of size segregation [36, 37, 40, 50]. Our experiments suggest that mixing of
particles by athermal self-assembly may indeed be feasible in a realistic granular system. In order to test the robustness of this self-assembly behavior, we repeated our experiments with steel beads (see SM), and found a qualitatively similar behavior. In particular, steel beads can also self-assemble into large domains of S1 crystal, but the overall crystallinity of the system typically remains lower.

### C. Liquid-crystal phase coexistence

In order to explore how the area fraction affects crystallization behavior, we perform a series of experiments at lower area fractions down to $\phi = 0.81$, while keeping the other parameters fixed. In Fig. 4, we show $q_4$ and $F_{S1}$ averaged over the last 50 minutes of each experiment as a function of $\phi$ (the total duration $t_{tot}$ varies for different experiments, see caption). We note a sharp increase of both observables between $\phi = 0.83$ and $\phi = 0.84$. The snapshots of the final configurations obtained in these cases (shown as insets) suggest that the system undergoes a liquid-solid-like phase transition as a function of the area fraction. For $\phi = 0.83$ the system exhibits an amorphous structure with no significant formation of S1 domains and the same holds also for lower $\phi$ (not shown). For $\phi = 0.84$, we found a single large domain of S1 particles coexisting with the amorphous phase. Finally, the last point at $\phi = 0.858$ coincides with the experiment discussed in the previous section, its final configuration is the one shown in Fig. 3e where the S1 domains are widespread all around the system. By looking at the final configurations in Fig. 4, it is possible to note a gradient of decreasing concentration of large particles from right to left (the same is also observed at $\phi < 0.83$). Of course, this does not happen in the S1 domain at $\phi = 0.84$ where the two species are perfectly mixed. We argue that this is a gravity-induced size segregation effect [37, 38, 41] due to the loss of horizontal calibration in the system. This was not observed in previous experiments done with the same setup [25], likely because they were performed at much lower vibration intensity. However, this experimental imperfection turns out to highlight the robustness of the S1 self-assembly process when the right geometrical conditions are matched. Indeed, S1 domains are able to homogenize the particle concentration even when gravity would promote size segregation. The systematic characterization of the competition between size segregation and self-assembly in granular binary mixtures goes beyond the scope of this paper and is left for future studies. On a practical level, this effect provides us with a qualitative criterion for deciding when to stop the experiments. If, after 20 hours of vibration, S1 self-assembly is not observed and relevant size segregation is found (with the consecutive formation of hexagonal clusters of small and/or large particles), we consider the experiment terminated. Otherwise, we let the experiment run until the formed S1 domains stop growing significantly.

The formation of a single crystal domain at $\phi = 0.84$ is strongly reminiscent of an equilibrium phase coexistence between a fluid and a crystal. This picture is further reinforced by the observation that the fluid phase in this system remains highly mobile, ruling out the possibility that the crystal has stopped growing because the system reached dynamical arrest. To demonstrate this, we compare in Fig. 5 the spatial variation in the averaged local $q_4$ with the one of large particle mobility for $\phi = 0.84$ and $\phi = 0.858$. Each map is performed by dividing the system into 15 mm side square subregions and averaging quantities over a 60-second-long high-rate acquisition (60 fps) done at the end of the experiments. To compute the local $q_4$, we simply average $q_4^m$ over each subregion, while as a measure of local mobility, we use the mean absolute displacement of the grains over a time interval of 0.17 seconds. From these maps, it is clear that, for $\phi = 0.858$, the mobility is quite homogeneous all over the system regardless of the local $q_4$. This means that the S1 phase is coexisting with other solid- or glassy domains (like the hexagonal one in the top-right corner of Fig. 3e). In contrast, for $\phi = 0.84$, the displacement map exhibits low mobility in the S1 phase (i.e. the region of large local $q_4$) and a much higher one in the surrounding fluid. The videos of the two experiments we provide in the SM further clarify the exposed picture. We point out that this phenomenology is consistent with the interpretation of the system at $\phi = 0.84$ forming a stable coexistence between a mobile fluid phase and a solid S1 cluster. We confirmed the reproducibility of this behavior in two additional runs (not shown) at $\phi = 0.84$, where we again observed the formation of low-mobile S1 clusters in coexistence a fluid.

In the remainder of the paper, we will further investigate the non-equilibrium phase coexistence between the granular binary crystal and the fluid phase found in our experiments by means of numerical simulations.

## III. NUMERICAL SIMULATIONS

As discussed in the previous sections, experiments to study the self-assembly of the S1 phase are extremely long and particularly hard to calibrate. In order to further explore the self-assembly process and better characterize the non-equilibrium phase coexistence between the granular liquid and the S1 phase, we now turn our attention to numerical simulations.

First, we study the numerical reproduction of the experimental setup using a realistic granular model that can be simulated through molecular dynamics based on the Discrete Element Method (DEM) [52, 53]. This analysis aims to show that our numerical simulations reproduce the S1 self-assembly process observed in the experiments.

Once the numerical model has been validated, we move to the implementation of the direct phase coexistence method [54–59] to systematically explore the behavior in the coexistence region. This will be done both with
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time averaged $q_4$ and $F_{S1}$ as a function of the area fraction for experiments with polyamide beads on a rough surface under sinusoidal vibrations. In all the cases $q = 0.5$, $x_S = 0.51$, $\Gamma = 1.79$ and $f = 53$ Hz. Time averages are performed over the last 50 minutes of the experiments. The total duration $t_{tot}$ varies for different experiments, we have $t_{tot} = \{20, 20, 20, 29, 89\}$ hours corresponding to the explored area fraction $\phi = \{0.81, 0.82, 0.83, 0.84, 0.858\}$. In the two insets, we provide the final snapshots for the experiments performed at $\phi = 0.83$ (top left) and $\phi = 0.84$ (bottom right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison between the maps of average $q_4$ and large particle displacement over 0.17 seconds for the experiments performed at $\phi = 0.858$ (panel a and b) and $\phi = 0.84$ (panel c and d). In the first case, the displacement field is homogeneous regardless of the local value of $q_4$, in the latter, low mobility regions correspond to the high $q_4$ cluster. The analysis is performed over a 60-second-long high-rate acquisition (60 fps) done at the end of the experiments. Each shown map results from an average of 3520 single-frame analyses (or double-frame in the case of the displacement).</figcaption></figure>

DEM simulations of the realistic granular model and with Event Driven Molecular Dynamics (EDMD) [60] of a simplified model which is able to capture the minimal ingredients of the granular dynamics under study [25, 61]. The direct phase coexistence method will allow us to identify the range of area fractions for which the liquid and the S1 phases are found in stable coexistence at different densities and, remarkably, at different granular temperatures.

## A. In silico replication of the experimental setup

DEM simulations consist of molecular dynamics that make use of realistic contact models to describe the motion of the grains during collisions with other grains or confining walls. This allows the reconstruction in silico of granular experimental setups, taking into account their real-world sizes and material properties. In our case, we use the Hertz-Mindlin model [53, 62, 63] for the contact dynamics as implemented in LAMMPS [62, 64, 65]. This model contains different parameters referring to the elastic and dissipative contribution of normal and tangential forces between colliding grains. The parameters related to elastic interactions are directly linked to the material properties (i.e. Young modulus and Poisson ratio) while the ones needed to model the dissipative terms must be fixed according to a calibration procedure. In our simulations, all the model parameters are tuned according to previous studies [65, 66] which revealed a good qual-
<figure><img src="image_7.png" /><figcaption>FIG. 6. DEM simulations of the quasi-2D vibrated setup with PBC. On the left, we show the initial and final configurations for $\phi = 0.83$ and $\phi = 0.85$ which evolve respectively towards a fluid and an S1 phase. On the right, we plot the BOOPs as a function of time for the same simulations. For all cases we use $N = 2000$, $q = 0.476$, $x_S = 0.5$, $\Gamma = 4.4$ and $f = 350$ Hz.</figcaption></figure>

itative agreement between experiments and simulations. We point out that we always use a Young's modulus lower than that of the real materials; this is a compromise that allows us to use a larger discretisation time step for the simulations thus reducing their computational cost. Our numerical setup mimics the experimental one: It is contained in a simulation box of height $h$ and width $L \gg h$ vertically confined by two horizontal plates. For the horizontal directions, we consider both Periodic Boundary Conditions (PBC) and hard walls. External energy is provided to the grains by imposing a sinusoidal vertical acceleration $a_z(t) = g\Gamma \sin(2\pi ft)$ to the boundaries of the system. In a preliminary analysis, we found that for the numerical system, a good driving condition to prevent grains from exploring the vertical direction when vibrated is $f = 350$ Hz and $\Gamma = 4.4$. It reasonably differs from the one used in the experiments since the simulated grains have material properties that do not perfectly match the real ones. Moreover, in this frequency range, we observed that the experimental system is subjected to the effect of vibrational modes of the bottom plate. Obviously, this does not occur in the simulation where hard boundaries are perfectly rigid. More details on the implementation of our DEM simulations can be found in the Supplementary Information of Ref. 25.

We initialize our system with a granular mixture having the desired geometrical parameters $\{q, x_s, \phi\}$ and ran-

dom horizontal velocities extracted from a Gaussian distribution. In Fig. 6, we show the initial and final configurations as well as the evolution of the BOOPs for two simulations with PBC at different area fractions $\phi = 83$ and $\phi = 0.85$. The system is prepared in an initial configuration where large and small grains are segregated into single-component hexagonal patches. For $\phi = 0.85$, once the vibration is turned on, the hexagonal domains melt and the two-grain species mix forming the granular S1 phase. Instead, at $\phi = 0.83$, the system evolves towards a homogeneous fluid phase. To get even closer to realistic conditions and test the robustness of the self-assembly process with respect to the system size, we now consider a larger system confined by hard walls at $\phi = 0.85$. In Fig. 7, we show the related numerical results for the final configuration, the BOOPs and the fraction of particles belonging to an S1 environment. First, we note that the formation of the crystalline phase occurs also in this case. Interestingly, the self-assembly of S1 in the simulations seems to be noticeably faster than in the experiment. If we compare the best experimental run in Fig. 3 with the two numerical ones at $\phi = 0.85$ in Fig. 6 and 7, we see that in the latter, the crystallization process is still ongoing at the end of the simulation even if a macroscopic fraction of the system is already in the S1 phase. We point out that in the simulations, hexagonal domains seem much easier to melt and this can depend on several
<figure><img src="image_8.png" /><figcaption>FIG. 7. DEM simulation of the quasi-2D vibrated setup with hard walls for a large system size *N* = 12000. a) Final configuration, b) time evolution of BOOPs and c) fraction of particles belonging to an S1 environment. The geometrical parameters are *φ* = 0.85, *q* = 0.476, *x*<sub>S</sub> = 0.5; the driving ones are Γ = 4.4 and *f* = 350 Hz.</figcaption></figure>

causes such as the fact that the simulation box is per-
fectly horizontal, the different driving parameters or the
specific granular interaction forces used in the numerical
model. The latter are realistic but still coarse-grained,
and hence cannot capture all the possible effects that are
present in the experiments.

To summarize, this numerical analysis confirms that
DEM modeling of granular materials can qualitatively
reproduce the S1 self-assembly found in the experimental
setup. This also means that this phenomenon is general,

rather than a consequence of some peculiarities of our
experimental setup.

B. Direct phase coexistence

In this section, we numerically explore the behavior of
the system in the coexistence region between the granular
liquid and the S1 crystal using the direct phase coexis-
tence method implemented both with DEM and EDMD
simulations.

In thermal systems, a common approach for deter-
mining equilibrium phase boundaries between coexist-
ing phases involves performing simulations in which the
two phases are in coexistence with each other. The abil-
ity of the two phases to exchange energy, volume, and
particles then naturally causes the two phases to evolve
towards thermal, mechanical, and chemical equilibrium.
This is the main idea behind the direct phase coexistence
method [54–57]. In this approach, one usually puts two
slabs of the different phases into contact in a periodic
simulation box elongated in one direction. This geomet-
rical setting usually stabilizes the interface between the
two phases if stable coexistence is reached. Depending
on whether the global value of the control parameter im-
posed on the system is inside or outside the coexistence
range, the system can either equilibrate towards a single
phase or towards stable coexistence. For coexistences in-
volving crystals, special care needs to be taken in order
to ensure that the crystal is not under any artificial strain
imposed by the boundaries of the simulation box [58, 59].
One crucial advantage of the direct coexistence method
is the fact that its implementation does not require an
assumption of thermodynamic equilibrium (in contrast
to e.g. free energy calculations [67]), and as a result, it
has also been applied to non-equilibrium systems such as
active matter [9, 68].

1. DEM Simulations

Our DEM numerical setup is the same as the one used
in Sec. III A but now the dimensions of the simulation
box are Lₓ × Lᵧ × *h* with Lₓ = λLᵧ where λ > 1.
The number N>L_ = N_S = 1800 and diameters σ_L = 2.5
mm, σ_S = 0.476σ_L of the grains are fixed while Lᵧ and
λ are chosen to satisfy the desired global area fraction
φ. Instead of initializing the system with a liquid and
a crystalline slab in contact, we prepare it with a single
crystalline block in the left-hand region of the box. The
initial lattice spacing *a*<sub>y</sub> of this S1 domain is such that it
occupies all the space in the *y*-direction. Then, depend-
ing on the global *φ* and the chosen area fraction of the
crystal slab, there will be a certain amount of free space
between the right edge of the initial crystalline block and
the one of the simulation box (see Fig. 8). Once ini-
tial velocities are randomly set according to a Gaussian
distribution and vertical vibrations are switched on, the
<figure><img src="image_9.png" /><figcaption>FIG. 8. Direct phase coexistence method for the DEM model. We show snapshots a) of the initial configuration, b) of a transient state, and c) of the non-equilibrium stationary state. We note that in the initial state, there is some free space between the right edge of the initial crystalline slab and the edge of the simulation box. In panel d), we plot the average $q_4$ as a function of time for five statistically independent runs. Simulations were performed with $N = 3600$, $\phi = 0.84$, $q = 0.476$, $x_S = 0.5$, $\Gamma = 4.4$ and $f = 350$ Hz.</figcaption></figure>

grains start to explore the available space. If $\phi$ belongs to the phase coexistence region, we expect from the equilibrium phenomenology that a part of the crystal will melt and the system will end up in a stable coexistence between a crystalline and a liquid slab. Despite our system being far from thermodynamic equilibrium, these expectations remain satisfied as we show for a specific case at $\phi = 0.84$ in Fig. 8. From the plot of the $q_4$ as a function

of time, we can see that, after a transient, the system reaches a non-equilibrium steady state with a stable degree of crystallinity. We remark that, in general, direct coexistence simulations lead to a deformed crystal with a different lattice spacing along $x$ and $y$ [58, 59]. To avoid this, one has to fine-tune the initial value of $a_y$ (or equivalently, the area fraction of the initial crystal slab) such that it is equal to the final $a_x$. All the results presented in this section have been obtained with $a_y = 1.0966\sigma_L$ corresponding to an initial crystal slab with area fraction 0.857. This gives the formation of an undeformed crystalline slab at coexistence (see the SM for our calibration procedure).

Since the system at $\phi = 0.84$ exhibits a phenomenology coherent with the one expected in the coexistence region, we now focus on the profiles of $q_4$ and area fraction around this state point. In Fig. 9, we plot the local BOOP $q_4(x)$ and local area fraction $\phi(x)$ along $x$ for $\phi = \{0.835, 0.84, 0.845\}$. For a given simulation, local observables are calculated based on bins of size $L_y \times \Delta x$ with $\Delta x = 2.2\sigma_L$ placed at $x \in [i\Delta x, (i+1)\Delta x]$ and averaged over 150 snapshots separated by 0.2 h in the steady state. The obtained profiles are then averaged again over three independent simulations. For each couple of profiles at a given global $\phi$, we observe a denser region with area fraction $\phi_{sol}$ corresponding to the S1 slab and a less dense one at $\phi_{liq}$ corresponding to the liquid phase. We note that the portion of the system occupied by the crystal increases with the global area fraction $\phi$ while $\phi_{sol}$ and $\phi_{loc}$ do not depend on it. This suggests the existence of an underlying lever rule for the density even if the system is not described by a free energy. Such behavior was also found and theoretically explained in active systems [69–71]. We can compute the average local area fraction in the two regions obtaining $\phi_{sol} = 0.855 \pm 0.001$ and $\phi_{liq} = 0.826 \pm 0.001$. To avoid the effect of the interface, $\phi_{sol}$ is calculated by averaging over bins with $q_4 > 0.96$ while $\phi_{liq}$ is calculated by averaging over bins with $q_4^i < 0.5$.

From this analysis, we conclude that the coexistence region of this granular liquid-solid-like phase transition is $\phi \in [0.826, 0.855]$. Remarkably, this is quantitatively coherent with the experimental results where we observed phase coexistence at $\phi = 0.84$.

## 2. EDMD Simulations

DEM simulations are particularly suitable to provide a direct link between real experimental setups and numerical simulations. Now we change perspective and approach the problem of granular phase coexistence with a model that provides a direct link with the reference equilibrium system where the self-assembly of the S1 phase was observed [45]. We consider the granular version of the non-additive hard-disk system described in Fig. 2a which is fully 2D and where collisions between particles are considered instantaneous events. The crucial differ-
<figure><img src="image_10.png" /><figcaption>FIG. 9. Horizontal profiles of $q_4(x)$ and $\phi(x)$ obtained with the direct phase coexistence method implemented for the DEM model. We report the results for three different imposed global area fractions $\phi = \{0.835, 0.84, 0.845\}$ and $N = 3600$, $q = 0.476$, $x_S = 0.5$, $\Gamma = 4.4$, $f = 350$ Hz.</figcaption></figure>

ences in the granular case are the coefficient of restituation $\alpha \in [0, 1]$ related to energy loss at collision and the explicit presence of a thermostat characterised by an inverse characteristic time $\gamma$ and a temperature $T_b$ (which, in order to be consistent with granular temperatures, we write in units of energy). This is implemented by imposing noise $\sqrt{2\gamma T_b/m_i \xi(t)}$ on the grains where $m_i$ is the mass of particle $i$ and $\xi = \{\xi_x, \xi_y\}$ represents a 2D delta-correlated continuous random variable distributed according to a Gaussian with zero mean and unitary variance: $\langle \xi(t) \rangle = 0$, $\langle \xi_i(t) \xi_j(t') \rangle = \delta_{ij}\delta(t-t')$. To summarize, in between collisions particle dynamics evolve following the Langevin dynamics [72]

$$ \dot{v}_i(t) = -\gamma v_i(t) + \sqrt{2\gamma T_b/m_i \xi(t)} \quad (2) $$

while particle collisions are governed by

$$ \begin{aligned} v'_{i} &= v_{i} + \mu_{ji}(1 + \alpha)(v_{ij} \cdot \hat{\sigma}_{ij})\hat{\sigma}_{ij} \\ v'_{j} &= v_{j} - \mu_{ij}(1 + \alpha)(v_{ij} \cdot \hat{\sigma}_{ij})\hat{\sigma}_{ij}, \end{aligned} \quad (3) $$

where $\mu_{ij} = m_i/(m_i + m_j)$ and $v'_i$ is the post-collisional velocities of particle $i$, while $\hat{\sigma}_{ij}$ and $v_{ij}$ are respectively the unit vector joining particles $i$ and $j$ and the relative velocity between them. For numerical simulations, we used a hybrid time-stepped/EDMD scheme as discussed in the SM. The model defined by Eqs. (2) and (3) has been extensively used to describe granular dynamics in different conditions [25, 73, 74].

We define the granular temperature of the system as its average kinetic energy

$$ T = \frac{1}{N_S + N_L} \frac{1}{2} \sum_{i=1}^{N} m_i \langle v_i^2 \rangle, \quad (4) $$

where $\langle \dots \rangle$ is an average over multiple uncorrelated configurations in the steady state. We also define per specie

granular temperatures $T_S$ and $T_L$

$$ \begin{aligned} T_S &= \frac{1}{N_S} \frac{m_S}{2} \sum_{i=1}^{N_S} \langle v_i^2 \rangle \\ T_L &= \frac{1}{N_L} \frac{m_L}{2} \sum_{i=1}^{N_L} \langle v_i^2 \rangle, \end{aligned} \quad (5) $$

where the sum runs only on either the small or large particles. We note that, for $\alpha = 1$, the system attains thermodynamic equilibrium at a granular temperature $T = T_b$ which is the same for the two disk species in virtue of the equipartition theorem: $T_b = T_S = T_L$. For $\alpha < 1$, the dynamics is out of equilibrium but can still reach a non-equilibrium steady state with an average kinetic energy that, in principle, depends on all the control parameters of the system $\phi, q, x_S, \alpha, \gamma, T_b$, as well as on the phase formed by the system. Out of equilibrium, energy equipartition between the two species is no longer guaranteed [75] and the velocity probability distribution function is not Maxwellian [73]. This model contains the essential ingredients of a granular system where the temperature $T_b$ accounts for the injection of external energy and $\alpha$ for the internal dissipation due to inelastic collisions. Here we are interested in exploring how the phenomenology of the liquid-solid-like phase transition under study is affected by the non-equilibrium mechanisms regulated by $T_b$ and $\alpha$. The parameter $\gamma$, as well as masses, diameters and compositions of the mixture, are kept fixed which allows us to define an energy scale $\epsilon = m_L(\sigma_L \gamma)^2$. Specifically, we choose $x_s = 0.5$ and $q = 0.476$ while the masses are fixed by requiring an equal mass density for both species.

Before exploring direct phase coexistence, we confirmed (not shown) that, for $\alpha = 1$, we recover the overall phenomenology observed in Ref. 45 where the equilibrium dynamics is performed in the microcanonical ensemble (i.e. without a thermostat) and that, for $\alpha < 1$, we can still observe the self-assembly of the S1 phase. For the direct phase coexistence simulations, we used the same general setup as the one described for DEM simulations but now our simulation box is fully 2D and we slightly changed the calibration procedure to set the correct initial lattice spacing to avoid deformed crystals at coexistence (see SM). In Fig. 10, we plot the $q_4(x)$ and the $\phi(x)$ profiles for a specific couple of $T_b$ and $\alpha$ observing results similar to the ones already found for DEM simulations thus validating this numerical scheme for further exploration of the granular phase coexistence. Then, we performed an extensive scan of direct phase coexistence simulations for $T_b \in [T_{min}, T_{max}]$ and $\alpha \in [0.85, 1]$. From each simulation, we extract $\phi_{liq}$ and $\phi_{sol}$ from the area fraction profiles so that we can build the phase diagrams reported in Fig. 11. We focus on the behavior obtained by varying $\alpha$ for fixed $T_b$ and vice-versa. In Fig. 11a, we can see that reducing $\alpha$ shifts the coexistence region towards larger area fractions thus promoting the stability of the granular liquid phase. Instead, in Fig. 11b,
<figure><img src="image_11.png" /><figcaption>FIG. 10. Horizontal profiles of $q_4(x)$ and $\phi(x)$ obtained with the direct phase coexistence method implemented for the EDMD model. We report the results for six different imposed global area fractions $\phi = \{0.85, 0.857, 0.863, 0.87, 0.877, 0.883\}$ and $N = 9600$, $\alpha = 0.94$ and $T/\varepsilon = 34.175$.</figcaption></figure>

we observe that reducing $T_b$ at fixed $\alpha$ has the opposite effect promoting the stability of the S1 phase with coexistence densities getting closer to the equilibrium values. It is interesting to note that two mechanisms that act in the same way from the energetic point of view (i.e. lowering the overall kinetic energy of the system) have a contrasting effect on the phase stability. In panels c) of the same figure we show the $T/T_b$ vs $\phi$ phase diagrams for all the explored values of $\alpha$ and $T_b$. This plot shows that the phase diagram in this representation remains qualitatively the same for different choices of $\alpha$. We can try to understand these observations by interpreting $T/T_b$ as a measure of the deviation from the equilibrium reference: at thermodynamic equilibrium, $T/T_b = 1$. Otherwise, the dissipative collisions ensure that $T/T_b < 1$, as the energy loss of the collisions is only partly compensated for by the thermostat. At the same time, we also know that dissipative collisions can act as an effective attraction potential [76]. This can explain the shift of the coexistence region as a function of $\alpha$ in Fig. 11a which is consistent with what is expected in an equilibrium system, where weak and very short-ranged attractions tend to shift the coexistence densities to higher densities [77]. Conversely, the results presented in Fig. 11b can be understood by considering the interplay between driving and dissipation mechanisms in the system. The Langevin bath drives the system towards equilibrium, while the dissipative collisions drive the system out of equilibrium. In the limit where $T_b \to 0$, the collision frequency between the particles approaches zero, whereas the characteristic timescale associated with the Langevin bath, $\gamma$, remains unchanged. Therefore, in this regime, the non-

equilibrium effects become negligible. On the other hand, at high driving temperatures $T_b$, the system is strongly driven out of equilibrium because the frequency of the dissipative collisions scales with $\sqrt{T}$.

To summarize, the phase diagram in Fig. 11c shows how the non-equilibrium mechanisms at play in the system change its phase behavior away from that of its equilibrium hard-disk analogue.

### C. Temperature difference between coexisting phases

Now that we have characterized the phase stability we focus on a remarkable non-equilibrium effect which can be observed both in the realistic and the coarse-grained model when the liquid and the S1 phase coexist. The usual picture at equilibrium is that there is no net flux of particles and heat between the two phases. This is because any macroscopic observation of a net current is intrinsically in contradiction with the time reversibility that characterizes equilibrium dynamics. From this, it is clear that the two phases must be in equilibrium with the same temperature and the same chemical potential. In out-of-equilibrium systems, this scenario can be violated: an open system where energy constantly flows in and out (as happens in granular or active matter) can sustain the presence of net currents while its statistical properties are stationary in time. This is a key difference between non-equilibrium steady states and equilibrium dynamics. As a result, the coexistence of phases with different average kinetic energies has been observed in active and granular monodisperse systems [9, 10, 16, 22, 27]. At the interface, this temperature difference can give rise to a net flow of energy from the hotter to the colder phase. Here we provide evidence of this distinctive out-of-equilibrium phenomenon in direct phase coexistence simulations of our granular binary system.

#### 1. DEM Simulations

In Fig. 12a-b we show the horizontal profiles of the granular temperature in the $xy$-plane for small and large grains obtained through DEM simulations of the same systems analyzed in Fig. 9. These quantities are defined as Eqs. 5 but considering only the horizontal projections of the grain velocities. First of all, we note the lack of energy equipartition: small and large grains have different average granular temperatures [75]. Comparing these profiles with the $q_4$ and area fraction ones (Fig. 9), it is clear that the crystal has a larger granular temperature than the fluid. We stress that this result presents a substantial difference with respect to previous studies of monodisperse systems [15, 16, 27], where the denser crystalline phase is usually “colder” than the coexisting less dense fluid one. In those examples of non-equilibrium phase coexistence, the origin of the cooling in the crys-
<figure><img src="image_12.png" /><figcaption>FIG. 11. a) Variation of the liquid (square) and solid (circle) coexistence density as $\alpha$ is varied in the EDMD model at $T_b/\epsilon = 25$. b) Variation of the liquid (square) and solid (circle) coexistence density as $T_b$ is varied at $\alpha = 0.94$. c) Variation of the liquid (square) and solid (circle) coexistence density as $T_b$ and $\alpha$ are varied. For each plot, the stars represent the measured equilibrium coexistence density values. $N = 9600$ in all simulations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Horizontal granular temperature profiles $T_{xy}(x)$ for small and large grains obtained with the direct coexistence method implemented with DEM simulations (same simulations of Fig. 9). These simulations have been performed in a square box with PBC for $N = 2000$, $q = 0.476$, $x_S = 0.5$, $\Gamma = 4.4$ and $f = 350$ Hz.</figcaption></figure>

talline phase was explained in terms of bistability of the grain-plate dynamics [15, 27, 78, 79] or by the common picture for which, in granular systems, denser regions present a higher rate of dissipative collision and in turn a lower average kinetic energy [16, 27].

The scenario observed in our DEM simulations needs a different explanation and we argue that it relies on the interplay between z-to-xy energy transfer and the local structure of the system. Indeed, apart from some specific mechanisms [80], it is reasonable to think that the main energy transfer from the plate to the granular system occurs in the vertical direction. Then, collisions between grains at different heights allow for a subsequent energy transfer from z to the horizontal plane. Such off-planar collisions are particularly favoured between spheres of different sizes because their centers are on average at different heights when they vibrate on a substrate. In turn, the S1 structure favours contact between small and large grains, so the z-to-xy energy transfer is expected to be more efficient in this phase than in the liquid. We believe that this is the most relevant effect explaining the coexistence between a “hot” crystal and a “cold” fluid in our quasi-2D granular system.

## 2. EDMD Simulations and Kinetic Theory

Coming back to the EDMD simulations, in Fig. 13 we plot the granular temperature profiles of small and big particles for the simulations analyzed in Fig. 10. The same phenomenology discussed for DEM simulations is recovered, namely a hot S1 crystal coexisting with a cold liquid. Additionally, here we note that the undercompressed crystal at $\phi = 0.883$ (i.e. at a area fraction below the melting point) is hotter than the coexisting crystals. This observation aligns with the usual picture for which, within a single phase, the granular temperature is mainly governed by the overall density [16, 27].

Despite the similar phenomenology obtained in EDMD and DEM, we cannot use the same argument about op-
<figure><img src="image_14.png" /><figcaption>FIG. 13. Horizontal profiles of $T_S$ and $T_L$ obtained with the direct phase coexistence method implemented with EDMD simulations. We report the results for six different imposed global area fractions $\phi \in \{0.85, 0.857, 0.863, 0.87, 0.877, 0.883\}$, $N = 9600$, $\alpha = 0.94$ and $T/\epsilon = 34.175$.</figcaption></figure>

timization of the energy transfer given by off-planar particle collisions since in the EDMD simulations, energy is directly injected in the $xy$ plane by the thermostat and collisions are dissipative events. Based on these premises, and knowing that the crystal is denser than the fluid, one would expect to find a crystalline phase colder than the coexisting fluid which is clearly in contradiction with our numerical results. To understand the behavior of the granular temperature observed in EDMD simulations we adapted the kinetic theory for granular mixtures developed by Barrat and Trizac [75] to our collisional model for non-additive disks defined by Eqs. (2) and (3) obtaining the following equations for the granular temperatures $T_i$ of each species $i \in \{S, L\}$:

$$ \frac{\sqrt{\pi\gamma(T_b - T_i)}}{2m_i} = \chi_{ii}\phi_i \frac{(1-\alpha^2)}{\sigma_i m_i^{3/2}} T_i^{3/2} + \chi_{ij}\phi_j \sqrt{\frac{\sigma_i}{\sigma_j^3}} \mu_{ji}^2 \left[ (1-\alpha^2) \left( \frac{2T_i}{m_i} + \frac{2T_j}{m_j} \right) + 4(1+\alpha) \frac{T_i - T_j}{m_j} \right] \sqrt{\left( \frac{2T_i}{m_i} + \frac{2T_j}{m_j} \right)}, \quad (6) $$

where $\phi_i = \pi\sigma_i^2/(4L_xL_y)$, and $\chi_{ij}$ are the pair distribution functions at contact between species $i$ and $j$. In this notation, when $i = S$, then $j = L$, and conversely. We point out that the latter quantities are not known a priori and generally depend on $T_i$, $T_j$, $\phi_i$ and $\phi_j$. We do not enter here into the details of the derivation, for which we refer to [75], but we only point out that the left-hand-side of the above equation accounts for the effect of the thermal bath on species $i$ while the right-hand-side expresses the average dissipation rate. Equating these two contributions gives the condition that must be fulfilled by the granular temperatures $T_S$ and $T_L$ in the non-equilibrium steady state. Of course, Eq. (6) does not directly provide us with an explicit analytical expression of $T_S$ and $T_L$ as a function of the control parameters. However, it tells us that the quantities that govern the granular temperatures of the two species in a binary mixture are the four products given by the terms $\chi_{ii}\phi_i$ and $\chi_{ij}\phi_j$. This already highlights the non-trivial interplay between local structure and global density in the energetic balance between forcing and dissipation. It is important to stress that Eq. (6) has been obtained following the standard scheme of kinetic theory [75, 81], then assuming molecular chaos and a Maxwellian velocity distribution for the grains. Both these assumptions are never exactly fulfilled in a granular material but, in dilute systems and for $\alpha$ close

to 1, the predictions of kinetic theory are usually in good agreement with simulation results [75]. Our simulations are performed at $\alpha = 0.94$ but, to probe the liquid-solid-like transition, they explore a range of very high densities. To test the applicability of these theoretical predictions to our simulations, we compared the value of the average granular temperature $T$ (Eq. (4)) obtained with kinetic theory (i.e. solving numerically the system of two equations given by Eq. (6)) with the one of EDMD simulations for a wide range of area fractions spanning from the dilute gas regime to the solid state. We point out that to find the solution from kinetic theory we must insert in Eq. (6) the values of $\chi_{ij}$ measured from simulations. The results of this analysis are reported in Fig. 14a and b where we first note a good agreement between kinetic theory and simulations for the whole range of explored area fractions. Then, we also observe that the granular temperature is always a decreasing function of the area fraction except around the liquid-solid phase transition. Here we see that there exists a range of area fractions where the crystal is hotter than the liquid despite being denser. Even if the kinetic theory cannot be used to directly predict the behavior of the system inside the coexistence region, we can take the values of the granular temperatures at the right and left edge of it as representative of respectively the liquid and the crystalline ones
<figure><img src="image_15.png" /><figcaption>FIG. 14. Kinetic theory for the total temperature in a binary mixture in the liquid and in the S1 crystal. a) Measured values of $\phi_i\chi_{ij}$ with respect to $\phi$ in the liquid and the solid. Inset is a zoom around the coexistence region. Note that in the solid $\chi_{ss} = 0$ because in a defect-free S1 configuration, pairs of small particles never touch. b) Comparison as a function of the area fraction, of the measured and predicted total temperature of the system (Eq. 6). The grey area marks the coexistence region between the liquid and the solid. Inset is a zoom around the coexistence region. c) Error between the measured temperature and the theoretical prediction $(T_{\text{measured}} - T_{\text{theory}})/T_{\text{theory}}$ as a function of the area fraction. Similar results are obtained for $T_S$ and $T_L$ alone. $N = 2 \times 223^2$, $\alpha = 0.94$ and $T/\epsilon = 34.175$.</figcaption></figure>

at coexistence. This means that kinetic theory is able to predict that the solid must be hotter than the liquid. An explanation for this behavior in the 2D model can be found by looking at the products $\chi_{ii}\phi_i$ and $\chi_{ij}\phi_j$ plotted in Fig 14c. Here we see that they are monotonously increasing with $\phi$ apart from around the coexistence region where the large-large and the large-small contributions experience a jump down and $\phi\chi_{SS}$ becomes zero (as we expect in the S1 structure). From the way these products enter the right-hand-side of Eq. (6) (which represents the dissipation rate), we can understand that the observed behavior of $\chi_{ii}\phi_i$ and $\chi_{ij}\phi_j$ suggests that the transition from the solid to the S1 structure can reduce the average dissipation in the system. Indeed, the dissipative contributions of LL- and LS-collisions are reduced while collisions between small grains are simply not possible in the S1 crystal.

The performed analysis suggests that within a single phase, the common scenario for which higher densities lead to an increase in the dissipation rate and in turn to a lower kinetic energy is applicable. However, when the fluid-crystal-like transition occurs, the effect of the local structure sets in and the granular temperature is primarily determined by the products between the global area fraction and the contact values of the pair distribution function. This allows the possibility of having denser granular crystals with a higher kinetic energy than less dense granular liquids.

Finally, to avoid possible confusion, we stress again here that our interpretation of phase coexistence at different granular temperatures in DEM and EDMD is based on two independent mechanisms that rely on the details of the system under study (e.g. off-planar collisions in DEM). There are many ways to introduce the effect of forcing, dissipation and confinement in granular models and we expect the temperature in the solid to be larger or smaller than the one in the liquid depending on these

choices. Nevertheless, we believe that the physical insights provided in this section offer two concrete examples of how the coupling between local structure and energy flows can originate a stable kinetic temperature difference between coexisting phases out of equilibrium.

## IV. CONCLUSIONS

In this study, we have shown experimentally that a granular binary mixture of spherical grains vibrated on a substrate undergoes a fluid-crystal-like phase transition between a disordered phase and a square binary crystal. The transition behaves remarkably similar to an equilibrium fluid-crystal transition, including the observation of a stable coexistence between the two phases. To further investigate this experimental observation, we used computer simulations of both a realistic granular model and a coarse-grained collision-based model. In both models, we again observe strongly equilibrium-like phase coexistences, obeying the lever rule. However, the similarity to equilibrium phenomenology is broken when looking at the granular temperature, which is different in the two coexisting phases. Using a kinetic theory for binary mixtures, we show how this temperature difference arises from the interplay between local structure and energy transfer mechanisms.

Our study demonstrates how crystal self-assembly can be controlled based on equilibrium predictions even in the presence of strong non-equilibrium effects. In particular, the case of binary granular crystals offers a way to realize mixtures with a spatially constant composition. Additionally, our results highlight how simple kinetic theories can shed light on the physics of non-equilibrium phase coexistences.
ACKNOWLEDGMENTS

We thank B. Darbois Texier for helping us find the
granular materials for the experiment and for the pre-
cious feedback on our early results, we also thank A.
Puglisi, A. Gnoli and E. Fayen for their invaluable sup-
port in setting up this project and S. Cabaret for the

design of the quasi-2D cell. This work has been car-
ried on with the support of Investissements d'Avenir of
LabEx PALM (Grant No. ANR-10-LABX-0039-PALM)
and funding from the Agence Nationale de la Recherche
(ANR), grant ANR-21-CE06-0039. Two of the authors
(RM and GF) partially wrote this paper while partici-
pating in a program at the Erwin Schrödinger Institute
(ESI) in Vienna.

[1] G. M. Whitesides and B. Grzybowski, Science **295**, 2418 (2002).

[2] A. Van Blaaderen, R. Ruel, and P. Wiltzius, Nature **385**, 321 (1997).

[3] Z. Cai, Z. Li, S. Ravaine, M. He, Y. Song, Y. Yin, H. Zheng, J. Teng, and A. Zhang, Chem. Soc. Rev. **50**, 5898 (2021).

[4] M. A. Boles, M. Engel, and D. V. Talapin, Chem. Rev. **116**, 11220 (2016).

[5] S. C. Glotzer and M. J. Solomon, Nat. Mater. **6**, 557 (2007).

[6] B. A. Grzybowski, K. Fitzner, J. Paczesny, and S. Granick, Chem. Soc. Rev. **46**, 5647 (2017).

[7] C. P. Royall, P. Charbonneau, M. Dijkstra, J. Russo, F. Smallenburg, T. Speck, and C. Valeriani, “Colloidal hard spheres: Triumphs, challenges and mysteries,” (2023), https://arxiv.org/abs/2305.02452:2305.02452 [cond-mat.soft].

[8] M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. **6**, 219 (2015).

[9] S. Mandal, B. Liebchen, and H. Löwen, Phys. Rev. Lett. **123**, 228001 (2019).

[10] L. Hecht, I. Dong, and B. Liebchen, Nat. Commun. **15**, 3206 (2024).

[11] A. K. Omar, K. Klymko, T. GrandPre, and P. L. Geissler, Phys. Rev. Lett. **126**, 188002 (2021).

[12] S. C. Takatori and J. F. Brady, Phys. Rev. E **91**, 032117 (2015).

[13] J. Bialké, J. T. Siebert, H. Löwen, and T. Speck, Phys. Rev. Lett. **115**, 098301 (2015).

[14] L. Caprini, E. Hernández-García, C. López, and U. Marini Bettolo Marconi, Sci. Rep. **9**, 16687 (2019).

[15] J. S. Olafsen and J. S. Urbach, Phys. Rev. Lett. **81**, 4369 (1998).

[16] A. Prevost, P. Melby, D. A. Egolf, and J. S. Urbach, Phys. Rev. E **70**, 050301 (2004).

[17] F. V. Reyes and J. S. Urbach, Phys. Rev. E **78**, 051301 (2008).

[18] A. E. Lobkovsky, F. V. Reyes, and J. Urbach, Eur. Phys. J. Special Topics **179**, 113 (2009).

[19] P. M. Reis, R. A. Ingale, and M. D. Shattuck, Phys. Rev. Lett. **96**, 258001 (2006).

[20] K. Roeller, J. P. Clewett, R. Bowley, S. Herminghaus, and M. R. Swift, Phys. Rev. Lett. **107**, 048002 (2011).

[21] J. P. Clewett, J. Wade, R. Bowley, S. Herminghaus, M. R. Swift, and M. G. Mazza, Sci. Rep. **6**, 28726 (2016).

[22] J. Olafsen and J. Urbach, Phys. Rev. Lett. **95**, 098002 (2005).

[23] G. Castillo, N. Mujica, and R. Soto, Phys. Rev. Lett. **109**, 095701 (2012).

[24] G. Castillo, N. Mujica, and R. Soto, Phys. Rev. E **91**, 012141 (2015).

[25] A. Plati, R. Maire, E. Fayen, F. Boulogne, F. Restagno, F. Smallenburg, and G. Foffi, Nat. Phys. **20**, 465 (2024).

[26] P. M. Reis, R. A. Ingale, and M. D. Shattuck, Phys. Rev. E **75**, 051311 (2007).

[27] P. Melby, F. V. Reyes, A. Prevost, R. Robertson, P. Kumar, D. A. Egolf, and J. S. Urbach, J. Phys. Condens. Matter **17**, S2689 (2005).

[28] Y. Komatsu and H. Tanaka, Phys. Rev. X **5**, 031025 (2015).

[29] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. **68**, 1259 (1996).

[30] T. M. Knowlton, G. Klinzing, W. Yang, and J. Carson, Chem. Eng. Prog. **90** (1994).

[31] P. Coussot, Rheometry of pastes, suspensions, and granular materials: applications in industry and environment (John Wiley & Sons, 2005).

[32] F. Li, P. Anzel, J. Yang, P. G. Kevrekidis, and C. Daraio, Nat. Commun. **5**, 5311 (2014).

[33] N. Boechler, G. Theocharis, and C. Daraio, Nat. Mater. **10**, 665 (2011).

[34] Q. Wu, C. Cui, T. Bertrand, M. D. Shattuck, and C. S. O'Hern, Phys. Rev. E **99**, 062901 (2019).

[35] A. N. Karuriya and F. Barthelat, Proc. Natl. Acad. Sci. **120**, e2215508120 (2023).

[36] S. Aumaître, C. A. Kruelle, and I. Rehberg, Phys. Rev. E **64**, 041305 (2001).

[37] J. C. Williams, Powder technology **15**, 245 (1976).

[38] K. Ahmad and I. Smalley, Powder Technol. **8**, 69 (1973).

[39] J. Bridgwater, Particuology **10**, 397 (2012).

[40] N. Rivas, P. Cordero, D. Risso, and R. Soto, New J. Phys. **13**, 055018 (2011).

[41] D. Huerta and J. Ruiz-Suárez, Phys. Rev. Lett. **92**, 114301 (2004).

[42] F. Nakai and K. Yoshii, “Reducing segregation in vibrated binary-sized granular mixtures by excessive small particle introduction,” (2024), arXiv:2405.07290 [cond-mat.soft].

[43] B. A. Grzybowski, A. Winkleman, J. A. Wiles, Y. Brumer, and G. M. Whitesides, Nat. Mater. **2**, 241 (2003).

[44] E. Fayen, A. Jagannathan, G. Foffi, and F. Smallenburg, J. Chem. Phys. **152** (2020).

[45] E. Fayen, M. Impéror-Clerc, L. Filion, G. Foffi, and F. Smallenburg, Soft Matter **19**, 2654 (2023).

[46] C. Likos and C. Henley, Philos. Mag. B **68**, 85 (1993).

[47] A. Kudrolli, Rep. Progr. Phys. **67**, 209 (2004).

[48] A. Baldassarri, A. Puglisi, and A. Sarracino, Comptes Rendus. Physique **16**, 291 (2015).
[49] I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. **78**, 641 (2006).

[50] P. Melby, A. Prevost, D. A. Egolf, and J. S. Urbach, Phys. Rev. E **76**, 051307 (2007).

[51] T. C. Hui, X. Zhang, D. Adiga, G. H. Miller, and W. D. Ristenpart, Lab Chip **24**, 966 (2024).

[52] P. A. Cundall and O. D. L. Strack, Géotechnique **29**, 47 (1979).

[53] T. Pöschel and T. Schwager, *Computational Granular Dynamics* (Springer, Berlin, 2005).

[54] A. Opitz, Phys. Lett. A **47**, 439 (1974).

[55] A. Ladd and L. Woodcock, Chem. Phys. Lett. **51**, 155 (1977).

[56] A. Ladd and L. Woodcock, Mol. Phys. **36**, 611 (1978).

[57] J. N. Cape and L. V. Woodcock, Chem. Phys. Lett. **59**, 271 (1978).

[58] J. R. Espinosa, E. Sanz, C. Valeriani, and C. Vega, J. Chem. Phys. **139** (2013).

[59] F. Smallenburg, G. Del Monte, M. de Jager, and L. Filion, J. Chem. Phys. **160** (2024).

[60] F. Smallenburg, Eur. Phys. J. E **45**, 22 (2022).

[61] R. Maire, A. Plati, M. Stockinger, E. Trizac, F. Smallenburg, and G. Foffi, Phys. Rev. Lett. **132**, 238202 (2024).

[62] https://docs.lammps.org/pair_granular.html.

[63] H. P. Zhang and H. A. Makse, Phys. Rev. E **72**, 011301 (2005).

[64] S. Plimpton, J. Comput. Phys. **117**, 1 (1995).

[65] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., Comput. Phys. Commun. **271**, 108171 (2022).

[66] A. Plati and A. Puglisi, Phys. Rev. E **102**, 012908 (2020).

[67] D. Frenkel and A. J. Ladd, J. Chem. Phys. **81**, 3188 (1984).

[68] B. Van Der Meer, V. Prymidis, M. Dijkstra, and L. Filion, J. Chem. Phys. **152** (2020).

[69] D. Evans and A. K. Omar, “Theory of nonequilibrium symmetry-breaking coexistence and active crystallization,” (2023).

[70] Y.-J. Chiu, D. Evans, and A. K. Omar, “Theory of nonequilibrium multicomponent coexistence,” (2024).

[71] S. Hermann and M. Schmidt, Phys. Rev. E **109**, L022601 (2024).

[72] A. Baldassarri, A. Barrat, G. D'anna, V. Loreto, P. Mayor, and A. Puglisi, J. Phys. Condens. Matter **17**, S2405 (2005).

[73] A. Puglisi, V. Loreto, U. M. B. Marconi, A. Petri, and A. Vulpiani, Phys. Rev. Lett. **81**, 3848 (1998).

[74] A. Sarracino, D. Villamaina, G. Costantini, and A. Puglisi, J. Stat. Mech. **04013**, P04013 (2010).

[75] A. Barrat and E. Trizac, Granul. Matter **4**, 57 (2002).

[76] R. A. Bordallo-Favela, A. Ramírez-Saíto, C. A. Pacheco-Molina, J. A. Perera-Burgos, Y. Nahmad-Molinari, and G. Pérez, Eur. Phys. J. E **28**, 395 (2009).

[77] L. Mederos and G. Navascues, J. Chem. Phys. **101**, 9841 (1994).

[78] W. Losert, D. G. W. Cooper, and J. P. Gollub, Phys. Rev. E **59**, 5855 (1999).

[79] J.-C. Géminard and C. Laroche, Phys. Rev. E **68**, 031305 (2003).

[80] These are: impacts with surface asperities (experiments) and coupling between rotational and translational dynamics mediated by the tangential friction with the plate (experiments and DEM simulations).

[81] N. V. Brilliantov and T. Pöschel, *Kinetic Theory of Granular Gases* (Oxford University Press, 2004).