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

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

► To cite this version:

A. Plati, R. Maire, F. Boulogne, F. Restagno, F. 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-05380354

HAL Id: hal-05380354

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

Submitted on 24 Nov 2025

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 Attribution 4.0 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: November 24, 2025)

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 non-equilibrium phase coexistences between a fluid and crystal at different granular temperatures [15, 27, 28], ii) a change from a two-step liquid-solid phase transition involving an intermediate hexatic phase to a one-step first-order-like scenario [29] 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 [30–36] 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 [37–42] 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 [43] 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 [44], 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, 45] and active [9, 10, 46] matter. Here,
equilibrium-like general principles able to predict their
occurrence and describe their characteristics are absent:
Each specific case needs to be analyzed to understand the
key physical ingredients that make the coexistence pos-
sible and determine its properties such as, for example,
which of the two phases is hotter. Extending a numer-
ical analysis we previously developed for monodisperse
systems [28], we show that in our systems the crystalline
phase has a higher granular temperature than the coex-
isting liquid one. We explain this observation by using
the relationship between the system's local structure and
energy transfer mechanisms. Remarkably, we find that
the kinetic 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

<figure><img src="image_2.png" /><figcaption>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 Lp = 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.</figcaption></figure>

$\rho = 1.14 \text{ g/cm}^3$. The system is confined in a quasi-2D
square box (height $h$ and width $L \gg h$) which is ver-
tically vibrated by an electrodynamic shaker (Brüel &
Kjær, LDS V400). Note that the particles are primar-
ily confined close to the bottom plate due to gravity,
rather than due to the presence of the top plate. The
main motivation for including the top plate is to pre-
vent rare fast-moving beads from escaping the system.
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 sandblasted (glass sand,
granulometry {425–600} µm) before anodization to im-
prove 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 projection 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 differ-
ently 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 genera-
tor (Keysight 33220A) is used to design a desired wave-
form $a_z(t)$ which fixes the vertical acceleration 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 mixture, we define the size ratio $q = \sigma_S/\sigma_L$, the fraction of small grains $x_S = N_S/(N_S + N_L)$ and the area fraction $\phi = (N_S\sigma_S^2 + N_L\sigma_L^2)\pi/4L^2$. Throughout all our experimental study we kept $q = 0.5$ and $x_s = 0.51$ fixed while varying $\phi \in [0.81, 0.858]$ which corresponds to a variation of $N_S + N_L$ in the range [2350, 2490]. Once a granular mixture 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 [47, 48]. 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 [47–49]. 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. To narrow our parameter space, we calibrated the driving parameters ($\{\Gamma, f\}$) by performing short self-assembly tests, optimizing for fast horizontal grain motion while minimizing small particle stacking. We found that this condition can be properly realized by setting $\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$ [48]. 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 gran-

<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 [47, 48]. 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>

ular matter [37–39, 41, 42, 50–53], 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 accounting 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
<figure><img src="image_4.png" /></figure>

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 full experimental 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.

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 average local square orientational order and a decrease in the hexagonal one. We point out that the main S1 formation is observed around the few largest clusters nucleated at relatively short times. The remaining time of the experiment 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 appearance 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 increase 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 formation 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 hexagonal 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 [38, 39, 42] and may therefore play a role in the formation of hexagonal domains. The loss of calibration also limits the duration of our experiments, making it difficult to observe the merging of all the different self-assembled S1 domains, which is expected to occur over long time scales. However, the important thing to note here is that these effects are rather marginal compared to the mixing of the two species within the multiple S1 domains which continue to grow despite the formation of competing structures. The reproducibility of this result
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 the mixing of particles of different sizes in a granular system is a paradigmatic problem with broad relevance across numerous applications, including pharmaceuticals, chemicals, food processing, mining, construction, and transportation, where controlling segregation and ensuring uniformity are essential for efficient and reliable performance [38, 40, 54–57]. As argued in a recent numerical study [43], the self-assembly of binary crystals in a polydispersed granular system can be a way to address the problem of size segregation [37, 38, 41, 53]. 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 [38, 39, 42] 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
<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>

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) [58, 59]. 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 [60–65] to systematically explore the behavior in the coexistence region. This will be done both with DEM simulations of the realistic granular model and with Event Driven Molecular Dynamics (EDMD) [66] of a simplified model which is able to capture the minimal ingredients of the granular dynamics under study [25, 67]. 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 mo-
<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>

tion 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 [59, 68, 69] for the contact dynamics as implemented in LAMMPS [68, 70, 71]. 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 [71, 72] which revealed a good qualitative 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 [25]. 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 bound-

aries 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 random 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 causes such as the fact that the simulation box is perfectly 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 coexistence method implemented both with DEM and EDMD simulations.

In thermal systems, a common approach for determining equilibrium phase boundaries between coexisting phases involves performing simulations in which the two phases are in coexistence with each other. The ability 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 [60–63]. In this approach, one usually puts two slabs of the different phases into contact in a periodic simulation box elongated in one direction. This geometrical setting usually stabilizes the interface between the two phases if stable coexistence is reached. Depending on whether the global value of the control parameter imposed 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 involving 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 [64, 65]. One crucial advantage of the direct coexistence method is the fact that its implementation does not require an

<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 $\phi = 0.85$, $q = 0.476$, $x_S = 0.5$; the driving ones are $\Gamma = 4.4$ and $f = 350$ Hz.</figcaption></figure>

assumption of thermodynamic equilibrium (in contrast to e.g. free energy calculations [73]), and as a result, it has also been applied to non-equilibrium systems such as active matter [9, 74].

### 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
<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>

box are $L_x \times L_y \times h$ with $L_x = \lambda L_y$ where $\lambda > 1$. The number $N_L = N_S = 1800$ and diameters $\sigma_L = 2.5$ mm, $\sigma_S = 0.476\sigma_L$ of the grains are fixed while $L_y$ and $\lambda$ are chosen to satisfy the desired global area fraction $\phi$. 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_y$ of this S1 domain is such that it

occupies all the space in the y-direction. Then, depending on the global $\phi$ 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 initial velocities are randomly set according to a Gaussian distribution and vertical vibrations are switched on, the 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$ [64, 65]. 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 [75–77]. 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
<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>

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 [48]. 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 differences in the granular case are the coefficient of restitution $\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 [78]

$$ \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} \mathbf{v}'_i &= \mathbf{v}_i + \mu_{ji}(1 + \alpha)(\mathbf{v}_{ij} \cdot \hat{\boldsymbol{\sigma}}_{ij})\hat{\boldsymbol{\sigma}}_{ij} \\ \mathbf{v}'_j &= \mathbf{v}_j - \mu_{ij}(1 + \alpha)(\mathbf{v}_{ij} \cdot \hat{\boldsymbol{\sigma}}_{ij})\hat{\boldsymbol{\sigma}}_{ij}, \end{aligned} \quad (3) $$

where $\mu_{ij} = m_i/(m_i + m_j)$ and $\mathbf{v}'_i$ is the post-collisional velocities of particle $i$, while $\hat{\boldsymbol{\sigma}}_{ij}$ and $\mathbf{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, 79, 80].

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 \mathbf{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 \mathbf{v}_i^2 \rangle \\ T_L &= \frac{1}{N_L} \frac{m_L}{2} \sum_{i=1}^{N_L} \langle \mathbf{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 model parameters $\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 [81] and the velocity probability distribution function is not Maxwellian [79]. 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 $\varepsilon = 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. 48 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
<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>

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, 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 [82]. 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 [83]. 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, 28]. 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.
<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).</figcaption></figure>

## 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 [81]. 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 crystalline phase was explained in terms of bistability of the grain-plate dynamics [15, 27, 84, 85] 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 [86], 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 (see e.g. [87] for a granular hydrodynamic description that takes into account this effect). 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
<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>

“hot” crystal and a “cold” fluid in our quasi-2D granular system. All of this mainly concerns the kinetic energy in the $xy$-plane. For the vertical direction, we verified that the temperature difference is barely visible, as expected given the dominant effect of plate-grain collisions on the $z$-component (not shown).

## 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, 28].

Despite the similar phenomenology obtained in EDMD and DEM, we cannot use the same argument about optimization 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. As we have shown previously for monodisperse systems [28], the behavior of granular temperature in coexisting systems can be understood by adapting via kinetic theory. In particular, we adapt the kinetic theory for granular mixtures developed by Barrat and Trizac [81] 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 [81], 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$. In a binary mixture, these terms are related to the frequency of collision between species $i$ and $j$. This already highlights the non-trivial interplay between local structure and global density in the energetic balance between forcing and dissipation. Given the few observed overlaps in the realistic system, $\phi_i\chi_{ij}$ likely act as proxies for collision frequency and energy dissipation, not only in EDMD but also in
DEM. However, it is not trivial to predict an average energy change per collision in the DEM model. It is important to stress that Eq. (6) has been obtained following the standard scheme of kinetic theory [81, 88], 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 [81]. 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, consistent with the observations in Ref. [28] on monodisperse systems. 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 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 primar-

ily 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 [28].

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 liquid-like 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 precious feedback on our early results, we also thank A. Puglisi, A. Gnoli and E. Fayen for their invaluable support in setting up this project and S. Cabaret for the
<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) 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. b) 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. c) 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. 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>

design of the quasi-2D cell. This work has been carried 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 participating in a program at the Erwin Schrödinger Institute (ESI) in Vienna. The authors acknowledge the use of the Ceres high performance computer cluster at the Laboratoire de Physique des Solides to carry out the research reported in this article.

[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, Reviews of Modern Physics **96**, 045003 (2024).

[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] R. Maire, A. Plati, F. Smallenburg, and G. Foffi, “Non-equilibrium coexistence between a fluid and a hot-ter or colder crystal of granular hard disks,” (2024), arXiv:2411.17531 [cond-mat.soft].

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[43] F. Nakai and K. Yoshii, Granular Matter 27, 7 (2025).

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

[45] L.-H. Luu, G. Castillo, N. Mujica, and R. Soto, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 87, 040202 (2013).

[46] L. Caprini, D. Breoni, A. Ldov, C. Scholz, and H. Löwen, Communications Physics 7, 343 (2024).

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

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

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

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

[51] A. Baldassarri, A. Puglisi, and A. Sarracino, Comptes Rendus. Physique 16, 291 (2015).

[52] I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. 78, 641 (2006).

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

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

[55] Y. Fan, K. V. Jacob, B. Freireich, and R. M. Lueptow, Powder technology 312, 67 (2017).

[56] D. Parsons, Powder Technology 13, 269 (1976).

[57] E. L. Paul, V. A. Atiemo-Obeng, S. M. Kresta, et al., Handbook of industrial mixing (Wiley Online Library, 2004).

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

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

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

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

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

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

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

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

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

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

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

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

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

[71] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolinteanu, 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).

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

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

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

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

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

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

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

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

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

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

[82] 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).

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

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

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

[86] 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).

[87] P. Maynar, M. García de Soria, and J. J. Brey, Physical Review E 99, 032903 (2019).

[88] N. V. Brilliantov and T. Pöschel, Kinetic Theory of Granular Gases (Oxford University Press, 2004).
# SUPPLEMENTARY MATERIAL

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

### DETAILS OF THE EXPERIMENTAL PROCEDURE

Here we provide additional details on the experimental setup. We recall that this work has been done with the same setup used in Ref. 1. The reader can also refer to its Methods section and Supplementary Information for details about image processing and acceleration measurements of the plate.

#### Horizontal calibration of the setup

The inclination of the plate is adjusted using three levelling feet that support the apparatus. Before each experimental session, we followed a calibration procedure consisting of different steps. First, a preliminary rough adjustment is made using an electronic level (Laserliner). Next, we vibrate a dilute monodisperse system of large grains (with diameters $\sigma_L$ = 2.5 mm for steel and $\sigma_L$ = 4 mm for polyamide) and analyze the time-averaged density field after a 20-minute data acquisition. The vibration parameters are set to $f$ = 120 Hz, $\Gamma$ = 1.48 for steel, and $f$ = 53 Hz, $\Gamma$ = 1.87 for polyamide, chosen to produce gas-like motion in the grains.

For a perfectly calibrated surface, we expect the density field to exhibit mirror symmetry along both the x- and y-axes, although uniformity is not expected due to dissipative effects near the walls. However, following the initial adjustment with the level, the density field consistently shows higher intensity in a particular direction, indicating the influence of gravity. To address this, we fine-tune the levelling feet based on the observed asymmetry and repeat the process to achieve a more symmetric density distribution. As mentioned in the main text, extended experiments may result in some loss of horizontal calibration, but this does not significantly impact the self-assembly process.

### EXPERIMENTAL CALIBRATION OF DRIVING PARAMETERS

In this section, we describe the initial experiments used to calibrate the driving conditions for the longer-term self-assembly runs. In these experiments, we explored the effect of the non-equilibrium parameters without substantially varying the geometrical constraints {$q, x_s, \phi$}. We considered two different materials for the grains and for each of them we varied the driving parameters to find good conditions for the self-assembly process. We used both polyamide (PROLABO) and steel (Marteau & Lemarié) beads. Note that in the main text, only the polyamide results are reported. For polyamide beads we have diameters $\sigma_S$ = 2 mm, $\sigma_L$ = 4 mm; the density is $\rho$ = 1.14 g/cm$^3$. For steel beads $\sigma_S$ = 1.2 mm, $\sigma_L$ = 2.5 mm, $\rho$ = 7.85 g/cm$^3$. and $L$ = 10 cm.

We first discuss the results of two different experiments, one with steel beads and one with polyamide beads; both of them were performed with sinusoidal vibrations. We recall here that the bottom surface used with polyamide grains has been sandblasted to increase its roughness. This is not possible with steel grains, which are harder and tend to smooth out the rough surface after repeated collisions, a problem that is particularly relevant in our case since we need to perform very long experiments. With steel beads we use a smooth surface as done in [1].

An important difference between the two considered cases is the range of driving parameters {$\Gamma, f$} we can impose on the system. As we found in a previous study [1], good driving conditions for granular self-assembly must be such that the system can effectively explore the spatial configurations in the horizontal plane while remaining relatively confined in the vertical direction to prevent small particles from jumping on top of each other and thus altering the effective 2D area fraction of the system. Hence, for both systems, we must find the combinations of $\Gamma$ and $f$ for which these requirements are satisfied. The key advantage of using a sandblasted bottom plate is that the presence of surface asperities significantly enhances the z-to-xy energy in plate-grain collisions [2, 3]. Thus, with polyamide beads on a rough surface, we are able to enhance the horizontal mobility using stronger vibrations that, in the case of steel beads on a smooth surface, would also imply a considerable exploration of the vertical direction in the grain dynamics. Quantitatively, for polyamide beads, we find an optimal driving condition for $\Gamma = \Gamma_p = 1.79$ at $f = f_p = 53$ Hz while for steel beads, we use $\Gamma = \Gamma_s = 1.28$ and $f = f_s = 120$ Hz. It is important to note that, although these driving conditions only differ by a factor $\Gamma_p/\Gamma_s = 1.40$ in acceleration, they correspond to displacement amplitudes of $A_p = g\Gamma_p/f_p^2 = 158\mu$m and $A_s = g\Gamma_s/f_s^2 = 22\mu$m for polyamide and steel particles respectively. This leads, in turn, to a change in the shaking energy of a factor $A_p^2 f_p^2/(A_s^2 f_s^2) = 10$. Finally, the values of $\sigma_S$ and $\sigma_L$ for polyamide
<figure><img src="image_16.png" /><figcaption>FIG. 1. a) Time evolution of $q_4$ and $q_6$ obtained with polyamide and steel beads under sinusoidal vibrations. For polyamide grains $q = 0.5$, $x_S = 0.5$, $\phi = 0.851$, $\Gamma = 1.79$, $f = 53$ Hz and the bottom surface is sandblasted. For steel grains $q = 0.48$, $x_S = 0.52$, $\phi = 0.855$, $\Gamma = 1.28$, $f = 120$ Hz and the bottom surface is smooth. b) Final configuration obtained with polyamide beads with several S1 patches widespread all over the system. c) Final configuration obtained with steel beads with large hexagonal clusters and a few small S1 domains.</figcaption></figure>

and steel beads are different but the imposed geometrical constraints {$q, x_S, \phi$} are approximately the same for the two experiments (see the caption of Fig. 1). In Fig. 1a-c, we plot the BOOPs (see the definition in the main text) as a function of time and provide a direct visualization of the final configurations of the system. The evolution of the BOOPs is qualitatively similar in the two experiments for the first 5 hours, exhibiting an increase of $q_4$ and a decrease of $q_6$. During this time, we observed the self-assembly of S1 domains in different regions of the system. After that, steel beads show an arrest of the S1 formation coinciding with the growth of hexagonal clusters of large particles which result in being particularly stable and hard to melt. For polyamide grains, we found instead a further increase of $q_4$ which is related to agglomeration and growth of the S1 patches. Here, no significant formation of hexagonal domains is observed. It can be concluded that the system with steel beads driven on the smooth surface at {$\Gamma_s, f_s$} is more prone to crystallization into hexagonal crystals than the one with polyamide beads driven on the rough surface at {$\Gamma_p, f_p$}. Here we try to explain this idea in terms of energy transfer mechanisms. We hypothesise that, in the system with steel beads, the energy is injected in the horizontal plane mainly by off-planar collisions and these are favoured when considering particles of different sizes. On the other hand, particles of identical size will have much less occasion of undergoing off-planar collisions with a relevant energy transfer between vertical and horizontal directions. Following this idea, it is reasonable to expect hexagonal clusters of grains of the same size to be particularly stable. As suggested by previous experimental studies [1], using granular mixtures with a higher fraction of small particles results in more homogeneous horizontal fluidization of the system. However, in our case, we cannot adopt this solution since we are interested in keeping $x_S \simeq 0.5$, i.e. close to the value that leads to S1 formation in the equilibrium reference system.

For the polyamide case, we can use a rough bottom surface, along with stronger vibrations, which strongly changes the situation. Indeed, grain-plate collisions alone are capable of transferring considerable energy from $z$ to $xy$, regardless of the local structure of the system. We conclude that a more efficient local injection of energy in the horizontal plane helps the system to remain homogeneously fluidized. Further experiments would be needed to definitively confirm this picture, but this goes beyond this study.

Experimentation with the shape of the driving signal showed that one way of improving the mobility of the steel beads was to use sawtooth driving signals instead of sinusoidal ones. Sawtooth vibrations are obtained imposing an acceleration $a_z(t) = g\Gamma[(ft \mod 1) - 1/2]$ so they are still parametrized by the frequency $f$ and the amplitude $\Gamma$. We found that with this type of driving it was possible to increase the strength of the vibrations without triggering a
<figure><img src="image_17.png" /><figcaption>FIG. 2. Results from two independent experimental runs with steel beads on a sandblasted surface under sawtooth vibrations. For both experiments $q = 0.48$, $x_s = 0.52$, $\phi = 0.856$, $\Gamma = 2.13$, $f = 110$ Hz. a) Time evolution of $q_4$ and $q_6$. b) Final configurations with S1 patches and no relevant presence of hexagonal clusters.</figcaption></figure>

relevant exploration of the vertical direction in the grain dynamics. This resulted in a stronger and more homogeneous horizontal fluidization of the system (but still weaker than the one obtained with polyamide beads on a rough surface). In Fig. 2, we show the results obtained with steel beads on a smooth surface under sawtooth vibrations with $\Gamma = 2.13$ and $f = 110$ Hz. We point out that in this case the S1 patches are more stable and the results are reproducible (we show data coming from two independent runs). However, when comparing the final configurations obtained with steel grains on a smooth surface and those obtained with polyamide grains on a rough surface, we always found larger and more widespread S1 domains in the latter case (compare for example Fig. 1b and Fig. 2b or 2c). Hence, in the main text, we only considered results obtained in the optimal experimental conditions we found for S1 self-assembly thus using polyamide beads on a rough surface under sinusoidal vibrations.

## CALIBRATION OF DIRECT PHASE COEXISTENCE SIMULATIONS

The direct coexistence method involves directly bringing two phases (performed as initial conditions) into contact around their respective coexistence densities and allowing the system to reach a steady state to measure its equilibrium coexistence properties. While this procedure is straightforward for fluid-fluid coexistence, additional care is required when at least one of the phases is a crystal. Specifically, it is important to ensure that the equilibrated crystal does not develop additional strain in the direction perpendicular to the interface due to boundary condition effects [4, 5]. Indeed, if the interface is perpendicular to the *x* axis, and the length of the simulation box in the *y*-direction is fixed, the crystal structure along the *y* axis cannot relax and notably, its lattice spacing will not change with time. Therefore, it is necessary to run multiple simulations with initial crystals of varying initial lattice spacing $a_y$, and identify the final steady system in which $a_y = a_x$. Here, $a_x$ is the measured lattice spacing, which can relax over time due to the fluid interface. This approach can be seen as a generalization of the stress-tensor-based method of Ref. 5. Although measurements of the crystal lattice spacing are somewhat more cumbersome (and likely less statistically accurate) than measuring the stress tensor, for non-equilibrium systems it has the advantage of avoiding the need for a robust method for measuring non-equilibrium stresses.

In DEM simulations, we start with a unique slab of S1 crystal and then let it melt until stable coexistence between the liquid and the crystalline slab is reached. We perform simulations with a fixed number of grains without varying
<figure><img src="image_18.png" /><figcaption>FIG. 3. Averaged $a_x/a_y$ as a function of $a_y$ for different area fractions in DEM simulations. Other fixed parameters are: $N = 3600$, $q = 0.476$, $x_S = 0.5$, $\Gamma = 4.4$ and $f = 350$ Hz.</figcaption></figure>

the number of columns and rows in the initial crystal. Then, we can tune $L_x$ and $L_y$ to satisfy the imposed global area fraction $\phi$. Here, the initial lattice spacing is fixed by $L_y$ since we start with an S1 domain that occupies all the space in the y-direction. Once the system has reached a steady state, we can measure the average ratio $a_x/a_y$ in the remaining crystal. We repeat this procedure for different $L_y$ and $\phi$ until finding the $L_y$ where $a_y/a_x = 1$, such that the crystal phase is undeformed. This process is repeated for different area fractions. In Fig. 3, we show for $a_y \sim 1.0965\sigma$ (corresponding to an imposed $L_y = 31.8\sigma$) we always satisfy this requirement having $a_y/a_x = 1$ with a maximum error lower than 0.004. We then kept $L_y = 31.8\sigma$ fixed for all the simulations discussed in the paper.

For EDMD, we use the following protocol. We run direct coexistence simulations for a sufficient duration to reach a steady state. A typical snapshot is shown in Fig. 4a. We directly measure the average lattice spacing in both the $x$ and $y$ directions, $a_x$ and $a_y$, in real space, as represented by the grey segments in Fig. 4a. The average lattice spacing can also be determined from the positions of the Bragg peaks in the structure factor. We verify that the final measured averaged $a_y$ is approximately equal to the value imposed initially (i.e. the number of layers or the orientation of the crystal did not change). We perform these measurements for various systems at the same global area fraction but with different initial $a_y$. This procedure results in Fig. 4c, showing the measured averaged $a_x$ as a function of $a_y$. The unstrained crystal corresponds to the crystal at which $a_y/a_x = 1$ which is found through a linear fit of $a_y/a_x$ as a function of $a_y$.

These simulations also allow us to obtain local area fractions, $\phi(x)$ through averaging multiple snapshots. This averaged profile is given in Fig. 4b. The average densities of the liquid and solid phases are determined by finding the mean values of the two plateaus in $\phi(x)$. Hence, for each simulation at a given $a_y$, the coexistence densities can be obtained, as seen in Fig. 4d. The correct ones correspond to those of the unstrained crystal, $a_y$ is such that $a_y/a_x = 1$ (stars), which was obtained in Fig. 4c. This method is then performed for different physical parameters to obtain the phase diagrams presented in the main text.

## HYBRID EVENT-DRIVEN/TIME-STEPPED MOLECULAR DYNAMICS

Since the potential of hard-disk particles is discontinuous, usual time-stepped molecular dynamics methods are not suitable for simulations of such systems. Instead, we use event-driven methods [6] where the time before the collision of two particles $i$ and $j$ denoted $t_{ij}^{col}$ can be analytically computed:

$$|\mathbf{r}_i(t_{ij}^{col}) + \mathbf{r}_j(t_{ij}^{col})| = \frac{\sigma_i + \sigma_j}{2}, \qquad (1)$$

with $\mathbf{r}_i(t)$ the position of particle $i$ at time $t$. For particles free flying or undergoing constant viscous drag, $t_{ij}^{col}$ can be found exactly from the initial velocities and positions of the particles. The viscous friction due to the thermostat is applied continuously on the velocities in between events. This exponential decay of the velocity leads to the following collision time between two collisions:

$$t_{ij}^{col} = -\log(1 - \gamma \delta t_{ij}) / \gamma, \qquad (2)$$
<figure><img src="image_19.png" /><figcaption>FIG. 4. Protocol used in EDMD simulations to obtain the solid and liquid coexistence densities. a) A typical coexistence configuration in its steady state. The segments in grey in the middle of the box show the numerically found values of the local lattice spacing (as shown in the Inset). N = 1944 b) The corresponding local area fraction obtained from an average of 150 (shifted) snapshots of the steady state pictured in panel a). c) Values of the average $a_y/a_x$ as a function of the initial (and because of stability, final) $a_y$. Each point corresponds to a different simulation with the same global area fraction but different initial $a_y$. The measurements are taken in the steady state and averaged over 10 configurations. The point at which $a_y/a_x = 1$ corresponds to an unstrained crystal. This informs us about the correct $a_y$ to impose at the beginning so that we can measure the correct coexistence densities. d) The measured liquid and solid densities (as seen in b)) for different initial $a_y$. The correct densities are the ones for which $a_y$ is equal to the one chosen from c). N = 9600, $q = 0.476$, $x_s = 0.5$ and $\alpha = 0.97$.</figcaption></figure>

with:

$$ \delta t_{ij} = \frac{-b - \sqrt{b^2 - v_{ij}^2 (r_{ij}^2 - (\sigma_i \sigma_j)^2)}}{v_{ij}^2} \qquad (3) $$

where $b = r_{ij} \cdot v_{ij}$ and $r_{ij}$ and $v_{ij}$ are respectively the relative position and velocity of particles i and j at the moment the subsequent collision time is computed.

The thermostat is modelled by a time-driven instantaneous change of the velocity of every particle each $\Delta t_{\text{noise}}$ (in between the particle-particle collisions) according to:

$$ v'_{i} = v_{i} + \sqrt{2\gamma T_{b} \Delta t_{\text{noise}} / m_{i}} \eta(t) \qquad (4) $$

with $\eta(t)$ a Gaussian white noise with zero average and unit variance. After the application of this noise, every collision time must be recomputed.

The fluctuation-dissipation theorem is respected (at equilibrium) for reasonable values of $\Delta t_{\text{noise}}$ because of the continuous nature of the damping [7]. We nonetheless chose $\Delta t_{\text{noise}} \ll \tau_f$, with $\tau_f$ the mean free time.

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

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

[3] G. Gradenigo, A. Sarracino, D. Villamaina, and A. Puglisi, Europhys. Lett. 96, 14004 (2011).

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

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

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

[7] L. Ma, X. Li, and C. Liu, Commun. Math. Sci. 15, 1171 (2017).