Observation of Liquid Glass in Molecular Dynamics Simulations

Molecular anisotropy plays an important role in the glass transition of a liquid. Recently, a novel bulk glass state has been discovered by optical microscopy experiments on suspensions of ellipsoidal colloids. 'Liquid glass' is a disordered analog of a nematic liquid crystal, in which rotation motion is hindered but particles diffuse freely. Global nematic order is suppressed as clusters of aligned particles intertwine. We perform Brownian dynamics simulations to test the structure and dynamics of a dense system of soft ellipsoidal particles. As seen in experiments and in accordance with predictions from mode coupling theory, on the time scale of our simulations rotation motion is frozen but translation motion persists in liquid glass. Analyses of the dynamic structure functions for translation and rotation corroborates the presence of two separate glass transitions for rotation and translation, respectively. Even though the equilibrium state should be a nematic, aligned structures remain small and orientational order rapidly decays with increasing size. Long-wavelength fluctuations are remnants of the isotropic-nematic transition.


I. INTRODUCTION
Colloids are ubiquitous in our everyday lives.Additionally, they play an important role as models, where they can be considered as 'big atoms' 1 .To date, most work exploring the structure and dynamics of colloidal suspensions has been performed with spherical particles.However, many physical and technological colloidal systems are non-spherical.Over the last two decades, the role of colloid form has come into the focus of scientific interest 2,3 .Shape matters especially at higher colloid concentrations, where the packing gets more intricate because of steric constraints.This is apparent already for ellipsoidal particles that are the simplest generalization of spherical bodies.For elliposidal colloids, theory and simulations predict a rich phase diagram of equilibrium fluid, crystal 4 and nematic phases 5 , plus additionally jammed [6][7][8] and glass states [9][10][11][12] due to the presence of structural and orientational correlations.The 3D structure of fluid dispersions of ellipsoids has been studied with confocal microscopy 13 .
Because the colloidal glass transition features many of the phenomena found in atomic and molecular systems, colloids have been used extensively as model systems to study the glass transition [14][15][16] .For a two dimensional system, glassy behavior in the orientation dynamics of anisotropic fluid was detected using event driven molecular dynamics simulations 17,18 .In three dimensions, decoupling between the translation and rotation degrees of freedoms was observed in molecular dynamics simulations of a binary fluid composed of dumbbell-shaped a) Electronic mail: mohammed-k-m.alhissi@uni-konstanz.deb) Electronic mail: andreas.zumbusch@uni-konstanz.de c) Electronic mail: matthias.fuchs@uni-konstanz.de particles where the rotation relaxation times were found to be smaller than the translation ones 19 .Using confocal microscopy of suspensions of prolate ellipsoidal colloids with an aspect ratio η = 3.5, we could recently show the importance of colloidal shape on the thermal, quiescent glass 20 .Liquid glass as a new glass state, predicted by theory 9 twenty years ago, was discovered.It is formed as intermediate material in an intriguing two-step vitrification process when increasing the volume fraction.Measurements of angular and translation time-dependent correlation functions demonstrated that when the packing density increases, the orientational motion freezes before the translation one.In between these two glass transitions, a liquid glass exists as a state where colloids can diffuse but do not de-correlate their local angular alignment.It was found that the global nematic order vanishes because ramified clusters of aligned particles interpenetrate and hinder each other 20 .
Glass states of ellipsoidal particles have been described in mode coupling theory (MCT) by Schilling and coworkers 9,[21][22][23] .MCT derives the structural arrest from the equilibrium static structure factors as input, and for larger aspect ratios predicted two different glass transitions for rotation and translation motion.The topology of state space was further elucidated by developing a schematic model of MCT 20 .It captures the generic coupling of two sets of degrees of freedom, translation and rotation ones, and captures their frequency-dependent interrelation in retarded friction kernels.The model predicts two generic MCT glass transitions, and a region in between that corresponds to the existence range of liquid glass.MCT makes specific predictions on the relaxation of the incipient glassy structure in the supercooled state 24 , that we will use to identify liquid glass in the present simulation work.
The formation of liquid glass is driven by the direc-arXiv:2401.01938v2[cond-mat.soft]14 May 2024 tional correlations accompanying nematic ordering.Yet, in some systems, these directional correlations preempt the ordering transition 20 , while in other systems, glass formation causes texture dynamics to arrest in nematic states.Highly elongated colloidal rods, for instance, show arrested space-dependent director fields in suspensions far in the nematic regime 25,26 .Static glassy structures obtained by driving particles in external fields [27][28][29][30] , and jammed states in athermal systems of ellipsoids 7,8 are other instructive examples for the importance of local alignment in amorphous solids made from anisotropic particles.A liquid glass-like state also exists in twodimensional colloidal films [31][32][33][34][35] , where it forms via the mutual steric hindrance of compact nematic domains.The nematic order of the two-dimensional domains extends over hundreds of particles, while it decays to zero fast in the three-dimensional liquid glass 20 .This suggests that ramified and fractal nematic precursors are the building blocks of liquid glass in bulk dispersions.Its mechanical properties should be intriguing as it is expected to flow like a liquid but to transmit torques like a solid even though it lacks nematic order.
Here, we report results of Brownian dynamics simulations shedding light on the microscopic structural correlations and the transport properties of liquid glass.The study is motivated by the intriguing absence of a clear observation of liquid glass formation in earlier simulations.To model the hard interaction potential of the ellipsoidal colloids found experimentally, event driven Brownian dynamics (EDBD) simulations 6,[36][37][38][39] were performed previously 10,20 .In these cases, only rather small systems of around 500 ellipsoidal particles could be simulated because of the numerical complexity to sample random motion without tolerating overlaps of anisotropic particles.In the two-dimensional simulations, an artificial dynamics using kinetic Monte Carlo (kMC) was employed for similar reasons 32 .The EDBD simulations reported the formation of a nematic phase as expected from equilibrium Monte Carlo studies.The formation of a nematic state could be suppressed by rough walls, which led to the formation of a (regular) glass 20 , viz. the state of arrested diffusion and rotation.The isotropic-nematic transition prevented the observation of a clear two-step relaxation in angular correlation functions as expected by MCT for liquid glass formation, and the identification of the two different glass transitions for rotation and translation required appreciable extrapolation of the rotation and translation diffusion coefficients to zero 10 .Similar extrapolations of the relaxation times of angular and translation functions was necessary in the kMC simulations in order to indicate the two glass transitions in the films 32 .No EDBD or kMC simulations have been reported in the liquid glass state, however.This has prevented simulations to test whether the dynamics in this state exhibits the decoupling of rotation and translation and follows the predictions of MCT.To overcome the limitations of the previous simulations, we first employ techniques that allow to simulate far larger ellipsoid systems than previously studied.Second, to be able to quench beyond the nematic state,we develop a soft repulsive potential that enables large steps in the effective volume fraction.It simplifies the established Gay-Berne potential 40 by only keeping the repulsion as has been done for spherical particles.For a test how close this approaches true hard interactions see e.g.Ref. 41 .A second aspect neglected in the previous EDBD simulations of ellipsoidal colloids 6,10,20,[36][37][38][39] concerns the anisotropic diffusion of the particles parallel and perpendicular to their main axis.In our work, we use the Brownian dynamics algorithm developed in Ref. 42 and take the rotation and translation diffusion coefficients from measurements of the liquid-glass forming colloids in dilute solution 20 .
The paper is organized as follows.In Sect.II the system as specified by its interaction potential is introduced.Sect.III summarizes the definitions of the structural and dynamic functions that are evaluated.In Sect.IV, aspects of the simulation methodology are given, and the central section, Sect.V gives the simulation results and compares them with experimental data where available.First, the isotropic and the nematic phases are characterized, and in Sect.V B the liquid glass is discussed.Sect.V C compares with MCT predictions and analyses the intermediate time dynamics of translation and angular correlation functions.The conclusions and a discussion of the results are given in Sect.VI.

II. MODEL
We study an isotropic dispersion whose constituents are soft ellipsoidal particles with a fixed aspect ratio.The interactions of the ellipsoids among themselves are governed by an anisotropic inverse potential.The study can build on approaches to molecular fluids, where efficient descriptions have been developed.In this section, the interaction potential is presented.
The interactions of anisotropic particles have been investigated for several decades 40,[43][44][45][46] .In 1981, Gay and Berne modified the overlap potential 43 introduced by Berne and Pechukas and presented their potential which has since been known as the "Gay-Berne potential" 40 .It is a generalization of the Lennard-Jones potential 47 in that it includes orientation dependent parameters.With it, the isotropic, nematic, and the other liquid-crystalline phases could successfully be predicted in molecular dynamics simulations 48,49 .However, the "Gay-Berne potential" is limited to anisotropic uniaxial particles of the same type.Berardi, Fava, and Zannoni extended it to describe the pair interaction of two identical biaxial molecules 44 and dissimilar biaxial molecules 45 .Their formulation describes the orientational degrees of freedom using the Euler angles.Berardi, Muccioli and Zannoni introduced another expression of the biaxial "Gay-Berne potential" using quaternions to describe the orientation dynamics of the anisotropic molecules 46 .This is more convenient for molecular dynamics simulations because the new potential with the quaternion formulation removes the singularity in the equations of motion resulting from use of the Euler angles to describe molecular rotations.
In this work we are mostly interested in an isotropic fluid at high densities.At such densities, the main forces contributing to structure formation are those due to the repulsive interactions of the particles' hard cores.In such a high-density fluid and since the system is far from the critical point, attractive interactions play negligible roles in stabilizing the fluid structure.Subsequently, we will therefore adopt a purely repulsive potential where the anisotropy enters through terms depending on the distance and the orientation of the ellipsoid.We are going to use only the repulsive part of the Gay-Berne potential developed in Ref. 46.The new potential will be referred to as a "Repulsive Gay-Berne potential (RGB)".Furthermore, we are going to study a one-component anisotropic fluid made up of identical uniaxial ellipsoids with an aspect ratio η = a/b.
The RGB potential describes the pairwise repulsive interaction of two ellipsoids, i and j, and is defined as: U (r, q i , q j ) = 4ϵ 0 ϵ(r, q i , q j ) σ 0 r − σ(r, q i , q j ) + σ 0 12 (1)  where ϵ 0 and σ 0 are an energy scale and a length scale respectively.r is a vector connecting the centers of masses of the two ellipsoids i and j.The 4-dimensional unit quaternion q i = (q i 0 , q i 1 , q i 2 , q i 3 ) specifies the orientation of the i th ellipsoid with the normalization condition |q i | = 1.The anisotropic contact term σ(r, q i , q j ) and the anisotropic interaction term ϵ(r, q i , q j ) are explained in detail in Ref. 46.In our work, these two terms are computed by setting the side-by-side, width-to-width, and end-to-end interactions to unity for a pair of ellipsoids with aspect ratio η = 3.5.
For the computation of orientation dependent quantities, transformations from the body frame into the spacefixed frame are required.The transformation defining the rotation for a unit vector from the i th ellipsoid body frame êb into the space fixed frame, where it becomes ês , is: where R is the rotation matrix written in terms of the quaternion components: Note that when considering the limit of isotropic molecules, Eq. 1 reduces to the repulsive part of the Lennard-Jones potential.

III. DEFINITION
In this section, we introduce the functions that describe the structure of our ellipsoidal fluid as well as the quantities that give us insight into the fluid's microscopic dynamics and correlations.

A. Statics
For the isotropic fluid, the structure of the fluid is identified by the radial distribution function g(r) and the isotropic structure factor S(q) defined respectively by where the brackets refer to the canonical average, r ij = r i − r j is the vector connecting the center-to-center distance between the two ellipsoids i and β, and the wave vector is q = (2π/L)n with L being the system length size and n a vector of integers.N is the total number of the ellipsoids in the fluid, and ρ = (N/L 3 ) is the particle number density.ρ(q) is the isotropic microscopic local density: The functions g(r) and S(q) do not give any information on the orientation correlations in the spatial and wave vector domains.Therefore, in addition to these two functions, we study two quantities that describe the spatial correlations of the orientation degrees of freedom.The orientation pair distribution function G n (r) is defined as where êi is a unit vector along the major axis of the i th ellipsoid expressed in a space fixed frame.P n is the Legendre polynomial of degree n.In the isotropic phase, G n (r) decays to zero while in the nematic and crystalline phases, it decays to nonzero values 48,50 .More precisely, for n = 2 this function decays to the quadratic equilibrium value of the scalar nematic order parameter for large distances 51 lim r→∞ G 2 (r) = S 2 .
Another interesting orientation-dependent quantity is defined in the wave vector domain.It is called the orientation structure factor 9 S lm (q): ρ lm (q, t) is the orientation-dependent microscopic local density which is defined in terms of the spherical harmonics Y lm (θ, ϕ) discussed in Ref. 52, While we only consider the density auto-correlations, viz.S lm , off-diagonal elements have also been measured 37 .We only compute the matrix element S 20 (q).This makes the azimuthal (the ϕ angle) dependent factor of S 20 disappear.The angle θ is the angle defined by the scalar product cos(θ) = q • ês where ês is a unit vector along the long ellipsoid diameter computed in the simulation box frame.After introducing S lm (q) and ρ lm (q), the equilibrium static structure of our fluid can fully be investigated.
An additional important quantity is the scalar nematic order parameter S which provides us with knowledge on the degree of the orientation order in our system.This scalar quantity is measured by computing the largest eigenvalue of the nematic order tensor 51 defined as: where u j α is the α component of the eigenvector ûj representing the orientation of the ellipsoid j in the simulation box frame.The largest eigenvalue then is where θ is the angle pointing away from the nematic director ⟨û⟩ = 1 N j ⟨û j ⟩.

B. Dynamics
In order to gain insight into the microscopic time evolution of the system, the dynamical quantities recording translation and orientation motions in the fluid need to be defined.These dynamical functions are defined in the temporal domain and tell us inter alia whether the system is in equilibrium so that we can compute the thermodynamic quantities.The first dynamical quantity is the well-known self intermediate scattering function 53 F s (q, t), This function probes single particle dynamics in the fluid and gives information on self correlation.In F s of Eq. (12) and in all other time-dependent correlation functions, time t ′ will be chosen large enough to ensure equilibration or aging will be discussed explicitly.No time-averaging over t ′ will be performed.In addition to self correlations, cross correlations among the ellipsoids in the fluid must be described.Therefore, we introduce the isotropic dynamic structure factor 53 S(q, t) , S(q, t) measures the collective dynamics and correlations of all particles in the fluid.To examine orientationdependent dynamics and temporal correlations, the following functions are adopted.The first one is the orientation dynamic structure factor 9 S lm (q, t) and the second function is the time-dependent orientation correlation function L n (t) defined in terms of the Legendre polynomial 54 with the same interpretation for the unit vector êj as in Eq. ( 7).The function S lm (q, t) measures the correlations of the coupled translation-rotation dynamics while the function L n (t) only measures the correlation of the rotation dynamics.As a special case for n = 2, the function L 2 (t) decays to the quadratic value of the equilibrium scalar nematic order parameter S for large timescales 51 , in that, lim t→∞ L 2 (t) = S 2 .Having introduced these functions, we can now probe the dynamics of the anisotropic fluid in detail.

IV. THE SIMULATION SETUP
We carry out Molecular Dynamics simulations using the Large-scale Atomic/Molecular Massively Parallel Simulator "LAMMPS" 55 .The Brownian dynamics of an ellipsoidal fluid is simulated at different densities using an algorithm implemented in LAMMPS.The integrator is known in LAMMPS, as the "fix Brownian/asphere", and is based on the work done by I. Ilie et al. 42 and S. Delong et al. 56 .Our system of interest is composed of 2197 ellipsoids with a fixed aspect ratio η = a/b = 3.5.In our Brownian simulations, the ellipsoids exist in a three dimensional box whose length size is L. Periodic boundary conditions are implemented in order to avoid any wallparticle interactions.The initial state of the system is prepared by uniformly distributing tiny isotropic particles on a cubic lattice.While the dynamics is integrated, the small particles are grown until they become ellipsoids with the required aspect ratio.After that, the equilibration phase was carried out taking into account the system density under study.For the isotropic states in section V A the system equilibration was recognized by the decay of the translation and orientation correlation functions to zero.For the nematic states, the equilibration was run until the scalar nematic order parameter only fluctuated around its average value.Once the equilibrations for the isotropic and nematic states were achieved, the measurement windows commenced directly for each state.For the liquid glass state the system was first equilibrated into an ergodic isotropic state at the small density Γ = 0.18.The liquid glass was then obtained by quenching the fluid into a high density Γ > 0.31.In order to relax the initial high forces, the system was then run until the translation correlations have relaxed.This gave the start of the measurement window for the liquid glass states in section V B. In the simulations, we set σ 0 = 2b = k B T /ϵ 0 = 1 where b is the ellipsoid short semiaxis, k B is the Boltzmann factor, and T is the temperature.The dimensionless time is t * = t/τ B where t is the time whereas the Brownian timescale is τ B = (2b) 2 /⟨D t ⟩ with ⟨D t ⟩ being the average translation diffusion coefficient.We refer to the dimensionless units with the star symbol.The timesteps used in the measurement phases are dt * = 10 −4 which is used in Sect.V A, and dt * = 5×10 −5 which is used in Sect.V B. The simulation time for Sect.V A is 10 6 dt * while the simulation time for Sect.V B is 60 × 10 6 dt * .The anisotropic translation and rotation diffusion tensors input in the simulations are the body-frame diffusion tensors measured in an experiment on dilute colloids made up of hard ellipsoids.All the experimental results or data in this research are taken from the experiment done by Roller et al. 20 .It is important to note that in the short time limit, the positioning accuracy of 20 nm in the xy-direction and 50 nm in the z-direction 57 and the limited statistics for long lag times contribute to the noise especially in the correlation data.The anisotropic translation diffusion tensor is D * t = diag(0.937,0.937, 1.125) while the anisotropic rotation diffusion tensor is D * r = diag(0.123,0.123, 0.476).This means that the average scalar translation diffusion is ⟨D * t ⟩ = 1.The dimensionless distance and wavenumber are r * = r/σ 0 and q * = qσ 0 , respectively.
As we use soft repulsive ellipsoids in the simulation, we adopt an effective density as parameter of the thermodynamic states.The density of the soft particles is notably temperature dependent.For the soft particles the effective density is defined as Γ = ρσ 3 0 (ϵ 0 /k B T ) 3/k .In our simulation, we set ϵ 0 = k B T and k = 12.Using the effective density is an approximation, as the potential in Eq. ( 1) has corrections to the isomorphs scaling at short distance 58 , yet Γ captures the dominant temperature dependence.With the above setup, the simulations run at different effective densities Γ.
Since most of the figures shown use reduced units, we are going to omit the stars from the dimensionless quan-tities for simplicity.When dimensional quantities are needed or discussed, their units will be stated explicitly next to them.Also, we will refer to the effective density only as the density.

V. RESULTS AND DISCUSSION
This section is divided into two parts: Identifying the equilibrium phases of the fluid (the isotropic and nematic) and finding the liquid glass as a nonequilibrium state in the soft ellipsoid fluid.For comparison with the experimental work by Roller, et al. 20 , experimental data are only available for the isotropic and the liquid glass states.

A. The Isotropic and Nematic Phases
In the simulations, soft ellipsoids with aspect ratio (a/b) = 3.5 are found to form isotropic states at densities Γ < 0.21.This can be inferred from the plots of the pair correlation function g(r) and the structure factor S(q) in Fig. 1 where there are no long range structures.In addition, the isotropic states manifest no nematic order which can be noticed from the decay of the orientation correlation in G 2 (r) to zero across the system length scales as noticed in Fig. 3. To inspect that there is no nematic order during the simulation time, the scalar nematic order parameter S is measured and shows a very weak and negligible nematic order (not shown).We attribute it to the finite simulation box size.In the plots of g(r) and S(q) we show comparisons between the soft ellipsoids at different densities Γ and the hard ellipsoids as experimentally measured at different volume fractions ϕ exp .We define a ratio of length-scales as mapping parameter λ = σ H 0 /∆ (where ∆ is some length scale) so that we can map the soft-ellipsoid structure onto the hard ellipsoid structure.Basically, we scale the x-axis of the soft-ellipsoid plots of the g(r) and S(q) functions by the factor λ until the structure of soft ellipsoids matches the structure of the hard ellipsoids.This means that for the structure plots of the hard ellipsoids we have 23 µm being the small diameter of the hard ellipsoid, since the hard ellipsoid plots are our reference plots.For the scaled soft ellipsoid plots, the best mapping parameter for the isotropic states is found to be λ = σ H 0 /σ S 0 = 0.83.Considering the translation degrees of freedoms, it is obvious from the g(r) and S(q) plots that the system structures in the simulation for the isotropic densities Γ = 0.16 and Γ = 0.18 are best mapped to the system structures in the experiment for the isotropic volume fractions ϕ = 0.40 and ϕ = 0.46, respectively.
At this point, the spatial orientation structure of the FIG.1: The translation structure functions g(r) and S(q) for two isotropic states in the simulation (solid and dashed lines) and their isotropic counterparts in the experiment (lines marked with squares and circles).ϕexp refers to the volume fractions in the experiment; one set of curves is shifted vertically for better visibility.For the experiment (hard ellipsoid) curve λ = 1.0 and for the scaled simulation (soft ellipsoid) curve λ = 0.83 (see the text).
isotropic states needs to be discussed.The upper panel in Fig. 2 presents the orientation radial distribution function G 2 (r) in an isotropic state.It is clear from the plot that the orientation correlations in the simulation, carried out with the anisotropic RGB potential, are more enhanced than the ones in the hard-ellipsoid experiment.Additionally, the G 2 (r) from EDBD simulations of hard ellipsoids 20 is included in Fig. 2. Its lower values reflect the fact that the soft ellipsoid neighbours align more strongly than the hard ones.This demonstrates the effect of the particle softness in the simulation which prevents producing unstable forces and torques at small length scales.The difference between the local angular structure in the soft-ellipsoid simulation and the hardellipsoid experiment results from the fact that we use an anisotropic potential in the simulation which enhances the local alignments of the soft ellipsoids.The orienta-FIG.2: The orientation structure functions G 2 (r) and S 20 (q) for isotropic states in the simulation (the solid and dashed lines) and their isotropic counterparts in the experiment (the lines marked with squares and circles).The G 2 (r) is shown for the isotropic state Γ = 0.18 and ϕexp = 0.46, respectively, and results from EDBD simulations of hard ellipsoids 20 are included as blue dash-dotted line.The S 20 (q) is shown for two volume fractions (one of them shifted vertically); see legends.The λ are the same as in Fig. 1.
tion structure factor S 20 (q) is shown in the lower panel of Fig. 2. The large wave vector oscillations show the discussed local differences, but the long-ranged orientation correlations become stronger on increasing the density equally in simulation and experiment.
When increasing the densities to larger than Γ ≥ 0.21, the soft-ellipsoid system transforms from the isotropic to the nematic phase.As there are no experimental data available for the nematic phase, we show the results of this phase only in the simulation.What characterizes the nematic states are the finite limit of the function G 2 (r) at far distances and its clear peaks as seen in Fig. 3. Additionally, the values of the nematic scalar order parameter S explicitly prove the existence of the nematic phase for Γ ≥ 0.21.They will be obtained from a finite size analysis performed in Fig. 16; see there for more details.It FIG.3: The simulation results for the spatial orientation correlation function G 2 (r) for densities given in the legend.It decays to zero and does not show a long range orientation order in the isotropic phase for Γ < 0.21.In the nematic phase, this function shows stronger orientation order that can be seen in the peaks at multiples of 2b.Also it does not decay to zero, but rather approaches the limit of S 2 for large distances; here S is taken from Fig. 16 (at particle number ≈ 275), see text for details.
is clear that when increasing the nematic density more and more, the degree of the nematic order increases as can be seen in the increasing values of G 2 (r ≫ 1) and S, that fulfill 51 G 2 (r → ∞) = S 2 .Probing the structures in the plots of the functions g(r) and S(q) in Fig. 4, we observe that the nematic states have more translation order than the isotropic states.The translation order is reinforced when the system density increases.The presence of a small peak in the S(q) plot around q ≈ 4.0 in Fig. 4 reflects the emergence of a weak translation order at that length scale.This weak translation order results from the strong orientation order at the short length scales in the orientation structure factor of the nematic states S 20 (q) around q ≈ 2.3 and q ≈ 3.8 in the same figure.Once again, the softness of the particles plays an important role in having such short and stable orders in the nematic states.An interesting feature of the function S 20 is that it informs about the long-wavelength orientation density fluctuations when looking at its small wave vector regime.It is clear that at small wave vectors in Fig. 4 the value of S 20 decreases when the fluid density increases.This is due to the fact that the amplitudes of the local orientation density fluctuations become smaller when the system moves away from the density at which the isotropic-nematic phase transition occurs (Γ ≈ 0.21).In the isotropic phase, the same phenomenon happens but in the opposite direction.When the system density approaches the isotropic-nematic transition, the value of S 20 at small wave vectors increases which is a sign that the magnitude of orientation density fluctuations increases so that these fluctuations dominate the whole simulation box at the phase transition.
Having discussed the structures of the isotropic and nematic states, we now discus their dynamics and show FIG.4: The pair correlation function g(r), the structure factor S(q), and the orientation structure factor S 20 (q) for different nematic states in the simulation; see legend.The x-axis is not scaled here.
the results of the translation and orientation correlation functions in both equilibrium phases.Figure 5 shows the self intermediate scattering function F s (q/λ, t) and the orientation correlation functions L 2 (t) and L 4 (t) in the isotropic states for the soft ellipsoids (simulation) and the hard ellipsoids (experiment) fluids.Notice that the time axis is in seconds, and that a single Brownian time τ B depending on the density via hydrodynamic interactions matches the dynamics of all time-dependent functions.The factor λ is present in the F s plots in order to have the correct mapping between the soft and hard ellipsoids at the respective length scales.It is obvious from Fig. 5 that the relaxation of the translation dynamics in the simulation agrees almost exactly with that in the experiment.This can be seen in the plots of the self intermediate scattering function F s for the two wave vectors (q/λ) = 6.42 and (q/λ) = 7.9.We repeat here that by setting λ = 1 we obtain the hard-ellipsoid value, while setting λ = 0.83 we obtain the scaled soft-ellipsoid value.The first vector (q/λ) = 6.42 corresponds to the peak in the isotropic structure factor and the second (q/λ) = 7.9 is a vector close to the peak seen in Fig. 1.Observing the decay of the self intermediate scattering function in the simulation and the experiment in Fig. 5, we deduce that the dynamics of the states Γ = 0.16 and Γ = 0.18 agree with the dynamics of the states ϕ exp = 0.40 and ϕ exp = 0.46 respectively.These agreements in the isotropic-phase dynamics between simulation and experiment have also been seen in the isotropic-phase structure functions g(r) and S(q) in Fig. 1.Also, the relaxations of the rotation degrees of freedoms depicted by the functions L 2 and L 4 in Fig. 5 show agreements and discrepancies in the rotation dynamics between the soft ellipsoids and the hard ellipsoids.In both systems, the soft ellipsoids and the hard ellipsoids, the rotation correlations seen in L 2 and L 4 tend to live longer and decay on rotation timescales larger than the translation timescales at which the F s functions decay.The rotation timescales are larger than the translation timescales by one to two decades in both the soft and hard ellipsoid systems.At the same time, the decay of the rotation dynamics in the hard ellipsoids sets in at earlier times than the one in the soft ellipsoids.An explanation for the discrepancy in the rotation correlation functions between simulation and experiment is the inevitable noise in the experiment that makes the correlation functions in the experiment decorrelate earlier than those of the simulation.Now, we concentrate our discussion only on the simulation results and show an overview of the dynamics in the isotropic and nematic states, as well as in the state of the liquid glass.Figure 6 depicts the plots of the translation correlation functions F s (q, t), S(q, t), and S 20 (q, t) at the wave vector q = 2.7, and the orientation correlation function L 2 (t).Notice the absence of the factor λ which means that the simulation plots are not rescaled because there are no experimental (hard ellipsoids) data.The wave vector q = 2.7 is close to the first small peak in the plot of S 20 (q) of the nematic states in Fig. 4. Figure 6 shows the isotropic states (Γ < 0.21 with solid lines), the nematic states (Γ ≥ 0.21 with dashed lines), and a liquid glass state (Γ = 0.33 with a solid line marked by circles).The fastest decay in the isotropic states corresponds to Γ = 0.16 while the slowest corresponds to Γ = 0.20.The fastest decay in the nematic states corresponds to Γ = 0.21 while the slowest corresponds to Γ = 0.25.In the supplementary materials, results are given for the structure and dynamics of the additional states Γ = 0.27, 0.29, 0.31.As expected, in each phase the relaxation timescale of the correlation functions increases when the density increases.At each density, the rotation timescale seen in L 2 (t) in Fig. 6 is remarkably larger than the translation timescales observed in F s (q, t), S(q, t), and S 20 (q, t) for wave vectors connected to the structural relaxation.The simulation time for the isotropic and nematic states stops once the states are equilibrated (at t = 100).We had to run much longer simulations for the state Γ = 0.33, and yet could not make angular correlations reach equilibrium.The simulation time for the state Γ = 0.33 ends at t = 3000 (see Fig. 6).In the nematic phase, the L 2 (t) function decays to a nonzero value which is the quadratic value of the equilibrium scalar nematic order parameter.Furthermore, the decays of the functions F s (q, t), S(q, t), and S 20 (q, t) in Fig. 6 for Γ = 0.21 are faster than the decays of the same functions at the isotropic state Γ = 0.20.This reflect the fact that when the particles cross the isotropic-nematic transition, their centers of masses can move more freely due to the orientation order which was absent in the isotropic state Γ = 0.20.However, this is only noticed close to the isotropic-nematic transition.At Γ = 0.33, a liquid glass state can be identified.From Fig. 6, the most obvious dynamical behavior of the liquid glass is the decay of the translation correlations in F s (q, t), S(q, t), and S 20 (q, t) around timescales t ≈ 700 while the orientation correlation in L 2 (t) persists with a very high amplitude.This reflects the fact that the particle rotation dynamics are arrested while the ellipsoid centers of masses are moving freely which enables the translation correlations to equilibrate.While in the nematic state, the arrest in the dynamic orientation correlations is in agreement with the measured nematic order, in the liquid glass, where S is very small (see Fig. 16), this relation clearly does not hold.
B. The Liquid Glass FIG.7: Structure functions sensitive to packing and alignment, g(r), S(q), G 2 (r), and S 20 (q), in the simulation at density Γ = 0.33 for three different simulations, i.e., three different waiting times t 1 w , t 2 w and t 3 w with t i w = 3000 × i where i ∈ 1, 2, 3.The decay of G 2 (r) to zero as r increases is a strong indication of the isotropic orientation correlations.
We now focus on the measurement phase of the liquid glass states at the density Γ = 0.33.The success in quenching into such a high density heavily relies on the nature of the particles used in the simulation.Introducing the particle softness by using the RGB potential helps to avoid the formations of unstable physical configurations upon quenching.Such softness does not generate any large forces that cease the MD simulation.
At the liquid glass densities, the anisotropic system becomes arrested in a metastable state identified by a plateau in the translation correlation functions ϕ t q (t) such as F s (q, t) and another plateau in the rotation correlation functions ϕ r n (t) such as L n (t).The possibility that at such high densities the system might develop long range translation and orientation orders at the same time, viz.enter the crystalline phase, can be inspected by looking at the structure functions in Fig. 7.The meanings of t 1 w , t 2 w and t 3 w in the figure refer to the waiting times for every plot to be produced in the simulation with t i w = 3000×i where i ∈ 1, 2, 3.This convention is adopted for Figs.7 to 10.Looking at the curves in Fig. 7, one can observe the lack of long range translation and orientation orders.In particular, G 2 (r) goes to zero as r increases which is a sign that the degree of alignment in the system at Γ = 0.33 is neither nematic nor crystalline 48,50 : Absence of crystallinity is obvious from the plots of g(r), G 2 (r), S(q) and S 20 (q) in Fig. 7 in the three different simulations (three different waiting times).An interesting feature at small wave vectors is a small peak in S(q) at q ≈ 2.7 which results from the emergence of the ori- entation correlations characterized by two distinct peaks in S 20 (q) at q ≈ 2.1 and q ≈ 4.2.This is reminiscent of the short range orientation order in the plot of S 20 (q) of the nematic phase in Fig. 4. The two peaks indicate that the development of the nematic phase is hindered by the slowing-down of the glassy dynamics.The presence of the two peaks has not been reported before, and it sheds light onto a strong local orientation order in the liquid glass at small wave vectors.We believe that the enhanced local angular structure results from anisotropic interactions.Contrary to the nematic states, the high density in the liquid glass state suppresses the free alignment of the ellipsoids and no substantial nematic order forms over the length scales in the simulation box.The scalar nematic order parameter S(t) for the three waiting times is shown in Fig. 8.It is obvious that there is no substantial nematic order in the three simulations.Probing the structure of the glassy state only does not enable us to identify the glass state.From Fig. 7 one cannot confidently distinguish between the structure of such a glassy state and the structure of the isotropic fluid discussed above.This aspect is typical of glassy systems.Hence the dynamics of the state Γ = 0.33 need to be discussed to obtain a conclusive explanation.
FIG. 9: Translation correlation functions in the simulation for Γ = 0.33 at q = 6.5.Fs(q, t) and S(q, t) are shifted vertically for a better view.The curves at three different waiting times t i w show that the translation dynamics do not manifest any aging; the translation degrees of freedom relax like in a fluid.
Concerning liquid glass dynamics, first the translation correlation functions shall be discussed.Fig. 9 shows the plots of the dynamic structure factor S(q, t), the self intermediate scattering function F s (q, t), and the dynamic orientation structure factor S 20 (q, t).These functions are computed for the wave vector of the main peak in the structure factor S(q = 6.5) as it can be seen in Fig. 7.It is clear that these functions do not show aging in the three simulations, viz.for the three waiting times.Interestingly, the translation dynamics illustrated in Fig. 9 behaves as the dynamics characterizing the supercooled hard-sphere fluid.After the initial diffusive regime, the dynamics relaxes slowly onto a fairly flat plateau and then decays below it.This functional form can be seen in the three translation correlation functions.The plateau is a clear signal of the existence of a cage.The final decay takes the form of a stretched-exponential curve with an initial power-law decay.It describes that the particles leave the cages which impede their translation motion for some intermediate time; see Sect.V C for more detail.
The orientation correlation functions exhibit a different type of dynamics.Figure 10 depicts how the orientation correlations evolve in the simulation.The figure shows L 20 (t).Because of their high plateaus, low orders of L n (t) would not give such clear indications for two step relaxation.As an example the plot of L 2 (t) is shown in the inset of Fig. 10.The reason for this is the tight arrest of the rotation motion of the ellipsoids.Hence, a higher order rotation correlation function (n = 20) is re-FIG.10: The orientation correlation functions L 20 (t) in the simulation at density Γ = 0.33 for three different waiting times t i w as labeled.Aging is observed by the dependence of the final decay on t i w .The aging seen in the three simulations evidences that the orientation motion can not relax in our simulation time window.The inset shows L 2 (t) which requires larger scale angular motion to decorrelate; thus aging arises in the first 2 % of the relaxation.
quired in order to resolve it.It is obvious from the figure that the ellipsoid rotation dynamics decay from the initial diffusive regime onto a plateau due to an angular cage which arrests the rotation motion.Later, the system dynamics fail to follow the α relaxation decay.The aging phenomenon causes the particle rotation dynamics to relax which is different from the case of the translation dynamics.The longer the simulation, the longer the plateau lives, which means that the particle orientations are arrested for longer times.The aging indicates that the system always finds a new metastable liquid glass state which has a deeper minimum in the free energy landscape of the system.While the orientation correlations are not able to relax, the translation dynamics has already relaxed much earlier around t ≈ 50 (see figures 9 and 10).The lifetime of the translation plateau is around two decades (and it is followed by a final relaxation), whereas the lifetime of the orientation plateau is around three decades.Hence, the observed type of cage is an orientation cage which allows the translation degrees of freedom to decorrelate but arrests the orientation ones.This confirms an assignment of the state as a liquid glass since the ellipsoids can move while their orientations remain arrested.The height of the plateau and the tightness of the cage will quantitatively be discussed in the framework of mode coupling theory for the glass transition in section V C.
Translation dynamics of the Brownian ellipsoids in the real space can be seen in the upper panel of Fig. 11, depicting mean squared displacements (MSD) for the isotropic, the nematic, and the liquid glass state.The bottom panel only shows MSDs for the liquid glass state Γ = 0.33 and its parallel (MSD ∥ ) and perpendicular The MSD for the liquid glass state Γ = 0.33 as well as its parallel (MSD ∥ ) and perpendicular (MSD ⊥ ) components.The plot of MSDq is obtained by using the self intermediate scattering function at q = 0.3 from the Gaussian approximation MSDq = (−6/q 2 ) ln(Fs(q = 0.3, t)).
(MSD ⊥ ) components.The reference point for the latter two is the ellipsoid's orientation at t=0.The bottom panel shows MSD q obtained from the self intermediate scattering functions F s (q, t) using the Gaussian approximation explained in the figure caption.The upper panel of Fig. 11 shows that the short time diffusion of all states is similar since the ellipsoids have not collided yet.Once the collisions start to take place around t = 2 × 10 −4 , the dynamics differs.At longer times arrest is seen as a weak plateau for the isotropic and nematic states (Γ = 0.16 − 0.25), and a clear plateau in the liquid glass state Γ = 0.33.At long times, particle motions again become diffusive for all states.The bottom panel of Fig. 11 shows that in the liquid glass state Γ = 0.33 the main contribution to the MSD plateau is due to the parallel motions of the ellipsoids.For small wave vectors, the self intermediate scattering function F s (q, t) can be used to obtain the mean squared displacement MSD q .The MSD and MSD q curves show similar behaviour over all time scales which confirms that at short and long timescales in liquid glass, the particle translation motions are diffusive, while at intermediate times 10 −1 < t < 10 caging arrests the translation motion of the ellipsoids.
FIG. 12: The structural functions, g(r) and S(q), at different liquid glass densities.λ is a scaling parameter which maps the simulation (soft-ellipsoid) curves to the experimental (hardellipsoid) curves.For the hard ellipsoids λ = 1.0 and for the soft ellipsoids λ = 1.04.The vertical dashed lines in S(q) mark the wave vectors used in Fig. 15, the vertical solid lines mark the q whose dynamics is analysed according to MCT in Fig.It is of interest to study the liquid glass at different densities.Figures 12 and 13 clarify the static structure functions of the liquid glass at different densities.In these figures the available experimental data are plotted.In the experiment, liquid glass was found at a volume fraction ϕ exp = 0.55.As before, the factor λ has the same meaning as in the isotropic phase in section V A, but here its best value for mapping the soft-ellipsoid liquid glass to the hard-ellipsoid liquid glass is λ = 1.04.From Fig. 12 it is obvious that the three isotropic states whose densities in simulation are Γ = 0.32, Γ = 0.33 and Γ = 0.34 map well onto the liquid glass state in experiment with volume fraction ϕ = 0.55.The λ value suggests that the FIG.13: The orientation structure functions, G 2 (r) and S 20 (q), at different liquid glass densities.λ is a scaling parameter which maps the simulation (soft-ellipsoid) curves to the experiment (hard-ellipsoid) curves.The upper panel includes the function G 2 (r) from a EDBD hard ellipsoid simulation.For the hard ellipsoid plots λ = 1.0 and for the soft ellipsoid plot λ = 1.04.In the lower panel, the vertical solid lines mark the q whose dynamics is analysed according to MCT in Fig. 18.
effective size of the soft and hard ellipsoids seems to be almost identical.The fact that both systems are at such high densities makes the packing constraints more important and similar, which results in λ ≈ 1.This agreement can be checked from the plots of the functions g(r) and S(q) in Fig. 12, which agree except for the small peak at q = 2.7 in the structure factor S(q) in the simulation.This small peak indicates that there is a weak short-range translation order at this length scale in the soft-ellipsoid fluid.This weak translation order in the simulation cannot be observed in the S(q) plot of the experiment (a weak indication of a shoulder can be noticed, however).It is directly related to the short range orientation order seen in the simulation plots of S 20 for the three densities; see Fig. 13.The two peaks in the simulation plot of S 20 (q) at q/λ ≤ 4.5 are a clear signal of the important role played by the ellipsoid softness discussed in the previous section.This short range orientation order is not seen in the experiment which reflects the fact that the orientation structure in the simulation does not exactly match the orientation structure in the experiment.A possible explanation for this is that the procedures of quenching the fluid into the liquid glass are not identical in simulation and experiment.In addition the ellipsoids in simulations interact via an anisotropic potential which clearly enhances the translation and orientation structures at short length scales.Another point to note in the simulation plot of S 20 is that its small wave vector limit decreases when the density increases.This is due to the fact that the isotropic-nematic transition density is further away at higher densities.Finally, the plot of G 2 (r) in Fig. 13 shows that the liquid glass state does not show any kind of global nematic order and that even the local nematic order becomes weaker as the density increases.This weakening of the global nematic order on increasing the density can also be observed in the global order parameter S shown in the lower panel of Fig. 8. Figure 13 includes the G 2 (r) of hard ellipsoids from the experiment and from the EDBD simulations of small systems 20 ; the latter show local alignment somewhat weaker but comparable to the soft ellipsoids, and otherwise are in a nematic state at the volume fraction where the dispersion forms a liquid glass.The G 2 (r) measured in the dispersion is very weak presumably for the reasons mentioned above when discussing S 20 (q).FIG.14: Dynamic structure and orientation structure factors, S(q, t) and S 20 (q, t), from the simulation at different densities in the liquid glass state for q = 6.5.
The dependence of liquid glass dynamics on the density is illustrated in figures 14 and 15. Figure 14 only contains simulation curves of the functions S(q, t) and S 20 (q, t) at q = 6.5 which is the unscaled location of the main peak in the simulation structure factor S(q).The dependence of the translation dynamics on the density only appears during the lifetime of the cage for every dynamic function (S(q, t) or S 20 (q, t)).Unexpectedly, the translation correlations increase when the density increases only at timescales starting at the β-relaxation regime until t = 50 in the α-relaxation regime.After that, the density dependence disappears and the translation dynamics for the three states Γ = 0.32, Γ = 0.33, and Γ = 0.34 behave identically, being characterized by a common timescale.This last finding does not agree with the MCT expectations since MCT predicts that the α relaxation timescale becomes longer when the density becomes higher 9,20,24 .However, some MCT predictions can still be seen, e.g. the increase of the plateau height as the density increases.In this context it is relevant that MCT considers the relaxation of an equilibrium fluctuation while the simulation records relaxations starting from a metastable state obtained after quenching.
FIG. 15: Self intermediate scattering function, Fs(q, t) (for q values marked in Fig. 12), and angular Legendre polynomial L 2 (t) and L 4 (t) at different densities in the simulation compared to the same functions in the experiment at volume fraction ϕexp = 0.55; parameters λ taken from Figs. 12 and 13.The time scale is mapped via τ B = 1.24 sec.
Figure 15 shows the comparison between simulation and experiment in terms of the translation dynamics illustrated by F s (q/λ = 6.42, t) and F s (q/λ = 7.9, t) and the rotation dynamics illustrated by the functions L 2 (t) and L 4 (t).The experiment could only probe the long time dynamics of the correlation functions in this case.The wave vector q/λ = 6.42 corresponds to the peak of the structure factor S(q) in Fig. 13 (the two wave vectors just mentioned are marked by the two vertical dashed lines in the S(q) plot).The behavior of the dy-namics of the function F s (q/λ, t) in the simulation is the same as the dynamical behavior of the functions S(q, t) and S 20 (q, t) when it comes to the density dependence and the MCT predictions.Figure 15 shows the best mapping in the dynamics when F s of the simulation is mapped at Γ = 0.33 onto the experiment at ϕ exp = 0.55.Both are liquid glass states.The Brownian timescale in this mapping is τ B = 1.24 sec which is much smaller than the Brownian timescales in the isotropic phase.The reason why the liquid glass Brownian timescale is much smaller than the Brownian timescales in the previous section (τ iso B = 140 sec and τ iso B = 220 sec in fluid states at ϕ exp = 0.40 and ϕ exp = 0.46, respectively) could be that at high densities the effective particle distances are smaller and thus the relevant energy scale in the softpotential is higher.Consequently, the corresponding collision rate would be higher and the particles would diffuse faster.
The close agreement observed in the translation dynamics is not seen in the orientation dynamics as the functions L 2 (t) and L 4 (t) illustrate.In the experiment, these two functions decorrelate to a lower value earlier than in the simulation.This is obvious from the decays of L 2 and L 4 in the experimental data around the timescale t = 50 while the decays of L 2 and L 4 in the simulation occur at much later timescales which cannot be observed in the time window of the two panels in Fig. 15 anymore.The reason for the lower correlation at longer times in the experiment presumably is the inevitable noise in the experimental environment which was also seen in the isotropic dynamics in Fig. 5. Importantly, however, also the experimental dynamic functions show the absence of a final relaxation of angular correlations.At all of the densities and in the translation and orientation dynamics of the liquid glass, we observe that the rotation correlation functions persist at least two decades longer than the translation one.
Finally, it is expedient to check how the nematic order changes as a function of the number of ellipsoids that form that nematic order.To this end, the simulation boxes are divided into a certain number of sub-boxes where each sub-box contains some number of ellipsoids.Then the nematic order parameter is computed in each sub-box and these values are averaged over the number of sub-boxes.The idea is to choose different numbers of sub-boxes leading to different numbers of particles per sub-box (or per cluster).Thereby, the average nematic order on different length scales is computed.Figure 16 clarifies how the magnitude of the nematic order changes as the number of particles per box increases (or as the number of sub-boxes decreases) for several densities.The figure shows that for all densities there is some varying degree of nematic order when the cluster contains a small number of particles.Once the particle number per cluster starts to increase, the nematic order in isotropic and liquid glass states starts to disappear quickly while the nematic order in the liquid crystalline states stays signifi-cant.This clearly proves the absence of compact nematic domains in liquid glass.Rather, confocal microscopy found ramified nematic precursors 20 .For the nematic states (at 0.21 ≤ Γ ≤ 0.25), there is a drop in S when it is calculated in the total simulation box.This indicates some domain structure which would require even longer simulation runs to homogenize throughout the system; thus the value of S for the eight-subbox system (average number of particles ≈ 275) is used when comparing to G 2 (r) and L 2 (t) in Figs. 3 and 6, respectively.Results of the cluster analyses for the states Γ = 0.27, 0.29, 0.31 are shown in the supplementary materials.
FIG. 16: Cluster analysis of the local nematic order in boxes of a given number of particles for different states (see legend) of the simulated soft ellipsoid system.Nematic order quickly vanishes for bigger clusters in isotropic and liquid glass states; see text for details.The value of S in the eight-box system (average number of particles ≈ 275), is used in Figs. 3 and 6.

C. MCT Analysis
Mode Coupling Theory (MCT) is the only microscopic theory to predict correctly many aspects of the structural relaxation close to the glass transition 9,24,59 .It is based on solving the equations of motions for the correlation functions describing the dynamics of molecular fluids.The main input of the theory is the static structure factor of the fluid under study.Here, we mention the parts of the theory most relevant to our work.The MCT-analysis we perform, studies the incipient frozenin correlations of the glass, how the system approaches it in the critical law (a power-law with exponent −a), and how it starts to decay during the von Schweidler law in the fluid (power-law with exponent b).In the glass, the long-time limit is a constant.The fitted MCT-laws are neither valid for short times, where colloids diffuse before colliding with others, nor for long times, where the liquid structure finally decorrelates following a stretched exponential.According to Franosch et al. 59 , the MCT FIG.17: The MCT fits to the translation correlation functions in the simulation at the density Γ = 0.33 for the different wave vectors q 1 = 3.2, q 2 = 6.5 and q 3 = 10.q 2 = 6.5 is the location of the maximum peak in the structure factor S(q).The S 20 (q, t) plot for the wave vector q 1 is shifted vertically for a better view.
predicts the time evolution of all correlation functions according to the following MCT solution where the positive sign refers to a glassy state, e.g. the rotation dynamics in liquid glass, while the negative sign refers to a fluid one, e.g. the translation dynamics in liquid glass.The label l becomes q (which refers to a wavenumber) when considering the translation dynamics.It becomes the Legendre-polynomial degree n when considering the rotation dynamics.The other quantities in Eq. ( 16) are the β-scaling function g ± (t/t σ ), which obeys 24,59 g ± (t/t σ → 0) → (t/t σ ) −a (17) for short and long rescaled times, respectively, and the correction functions These functions and parameters depend on the specific glass transition considered and hold for all correlators Φ l 59 .The amplitudes f l , h l , K l and R l are l-dependent.f c l gives an approximate estimate for the nonergodicity parameter at the singularity point (the glass transition density) while h l estimates the tightness of the cage arresting the dynamics.As we have two types of dynamics, we assert that there are two singularity points.The first is the rotation glass-transition density (or the fluid to liquid-glass transition density) labeled as Γ r c .The second singularity point is the translation glass-transition density labeled as Γ t c ; it leads from liquid glass to (regular) glass 20 .According to MCT discussion in Ref. 9 , using the density as a control parameter, the liquid glass appears before the translation glass (Γ r c < Γ t c ).We call all other parameters in the MCT solution, Eq. ( 16), the fitting parameters.The separation fitting parameter is defined as σ = σ 0 (Γ − Γ c )/Γ c , where Γ c is a singularity point at which the MCT solutions bifurcate (σ = 0).σ 0 is some constant, and the scaling time is given by t σ = t 0 /σ (1/2a) where t 0 is a system-specific cross-over time.The remaining parameters depend on the correlation function, and their values will be mentioned below.We implement the MCT analysis at the liquid glass density Γ = 0.33 since it exhibits a time-dependence close to the experiments.Figure 17 shows the results of the best fits of the translation dynamics to the translation part of Eq. ( 16) at Γ = 0.33. Figure 18 shows the best fits of the orientation dynamics to the orientation part of Eq. ( 16) at the same density.
For analyzing the translation correlation functions, the fitting parameters are set equal to the ones for the supercooled hard-sphere fluid 59 except for the exponents (a    The exponents a and b are universal, meaning that they are the same at all length scales while fitting the translation correlation functions (Fig. 17).They are different for the orientation correlation functions but are constant at all Legendre orders n while fitting the orientation correlation functions L n (t) (Fig 18).From Fig 17 we can notice that the translation correlation functions show agreement with MCT predictions around the turning point at the end of the plateau region (t ≈ 1) until some late timescales in the α relaxation regime, depending on the length scale we observe.Specifically, the MCT solution shows more agreement with the softellipsoid translation dynamics at small wave vectors (note the curves at q 1 = 3.2 and q 2 = 6.5 in Fig. 17) than the translation dynamics at the large wave vectors (note the curve at q 3 = 10.0 in the same figure).The discrepancy between theory and simulations at the large wave vectors lies in the β relaxation process.It is noteworthy that the long time decay in the translation dynamics is slower than predicted by MCT.It is expected that the discrepancies become smaller when the translation glass transition is approached, i.e. when moving closer to the translation bifurcation point of the MCT labeled as Γ t c .The values of the nonergodicity parameters f c q and the critical amplitudes h q at different length scales are shown in the left panel of Fig. 19.The nonergodicity parameter gives an estimate of the frozen-in amplitude in the respective translation correlation function at the translation-glass singularity Γ t c .The decrease of f c q values with the increase of the wave vectors indicates that the translation glass structure allows more motion at smaller length scale.The critical amplitudes h q for the three translation correlation functions are almost identical.This means that the tightness of the cage estimated by the three translation function at a certain wave vector is almost equivalent.The h q values increase as the wave vectors increase indicating that the glass structure of the ellipsoids gets arrested by a rather tight cage.Neglecting quantitative differences and because most of our fitting parameters were taken from the ones of the supercooled hard-sphere fluid, we find that the translation dynamics of the soft-ellipsoid liquid glass state behaves similarly to the translation dynamics of the supercooled hard-sphere fluid.
The fits of the orientation dynamics depicted by L n (t) in Fig. 18 are discussed for two different orders n = 10 and n = 20.The higher the order, the weaker the orientation correlation becomes with time (even in liquid glass).It is clear from the figure that the MCT solution agrees with the dynamics in the β-relaxation regime.The solutions do not show any agreement with the curves in the region where the aging-relaxation starts (t ≈ 100).The fact that at a comparable time for all L n (t) no agreement between the MCT solution and the L n curves can be found any more, is a signal of the break-down of the theory.The simulations clearly showed the aging of the orientation dynamics; see Fig. 10.Here, we analytically assert that the cooperative rotations within their orientational cage in the liquid glass state at Γ = 0.33 move according to the MCT intermediate-time dynamics of glass.
The nonergodicity parameter f c n and the critical amplitude h n are shown in the right panel of Fig. 19.The value of f c n gives an approximate value for the parameter at the rotation singularity point Γ r c .The f c n values become smaller when the order n increases.This tells that the orientation correlation becomes weaker when the rotation correlations are probed more locally.The critical amplitude h n illustrates how much the cage affects the orientation dynamics.At higher order n, the orientation cage of the ellipsoids exhibits larger motional amplitudes.

VI. CONCLUSIONS
We performed large scale Brownian dynamics simulations of ellipsoidal particles (at aspect ratio 3.5) including dense systems quenched into a partially arrested isotropic state.It is the first realization in simulations of the recently discovered liquid glass state.Based on our timeand wave vector-windows that appreciably extend the experimentally accessible ones, we tested the predicted decoupling of translation and orientation motion.Glassy arrest of the angular correlations was evidenced observing aging, while translation correlations relaxed in (supercooled) equilibrium and at least a factor of one hundred faster.Our simulation results compare well to the experimentally measured structure and to the translation motion in the colloidal dispersions where liquid glass was first observed 20 .While the long-wavelength static orientational correlations agree between simulation and experiment, the local alignment of the elongated particles showed differences.Presumably as a consequence, the time-dependent orientational correlations agreed in their final relaxation or when arresting, but showed differences for shorter times.In the simulations, angular correlations remained higher for longer, which might be caused by the fact that the pairwise interaction in simulations is anisotropic.
By a detailed analysis of the power-law relaxation in the intermediate time dynamics, we tested whether mode coupling theory correctly describes the diffusion and orientation dynamics of the liquid glass state.Overall we found good agreement of the time-dependent structural functions over typically four orders in time-variation.Note, that we ignored two failures of the microscopic MCT.First, it is well-known that the MCT glass transitions are rounded 24 .We identified the liquid glass state by aging thereby bypassing the question whether the final relaxation time is actually infinite or just longer than the simulation time.The second aspect concerns the more specific prediction in Ref. 9 that translation motion in liquid glass exhibits a small non-ergodic amplitude, while we find it to vanish.We follow the interpretation of the authors and previous tests in two [31][32][33][34][35] and three 20 dimensions to identify liquid glass with the state where rotation motion freezes while the structure on the average particle separation remains fluid.In the orientational dynamics we could only test the decay onto an intermediate time plateau because the final aging-induced relaxation is beyond the theory.Quenching into metastable states did not allow us to explore the change of the dynamics inside the liquid glass states as predicted by theory.Quenching to final glass states at higher effective density did not slow down translation motion.Presumably, the system falls out off-equilibrium during the quenching process in a way not affected by the final density aimed for.
The strong orientational fluctuations that we recorded in the static angular structure factor support the prediction by Letz, Latz and Schilling 9 , that liquid glass originates in nematic fluctuations, which become the entities of the mutual hindrance and of glassy arrest.During liquid glass formation, the local alignment of ellipsoids gets frustrated on intermediate length scales and gives ramified structures, viz.differently oriented clusters.The emerging nematic precursors contain more and more particles, and it is the mutual obstruction of such cooperative regions that leads to the formation of liquid glass.For future work it would be interesting to study larger simulation systems and and include hydrodynamic interactions in order to slow down the formation of global nematic order more so that the transition to the liquid glass state, where diffusion proceeds but angular motion freezes, can be studied using computer simulations.

VII. SUPPLEMENTARY MATERIALS
Refer to the supplementary materials which show the comparisons of the structure and dynamics in the liquid glass and the nematic states preceding the nematic-toliquid-glass transition.

Supplementary Materials
Observation of Liquid Glass in Molecular Dynamics Simulations Mohammed Alhissi, Andreas Zumbusch, Matthias Fuchs Results for the densities Γ = 0.27, 0.29, 0.31 In the following, we report results for the states with Γ = 0.27, 0.29, 0.31.At Γ = 0.27 and Γ = 0.29, the system is in a nematic state, whereas the state at Γ = 0.31 is a liquid glass state.Structure and dynamical functions are those defined in section III of the main article.The same simulation setup as described in section IV of the main article was used, except for the state with Γ = 0.29.The latter required longer equilibration and measurement times than the other nematic states.Its equilibration ran until the scalar nematic order parameter fluctuated around its average value.The measurement time of the state Γ = 0.29 was 10 × 10 6 dt * with dt * = 0.0001.Since none of the plots in these supplementary materials contains experimental data, the mapping parameter λ does not appear in the figures below.
The pair correlation g(r) and the structure factor S(q) functions giving insight into the translation order in the soft ellipsoid fluid for the states Γ = 0.27, 0.29, 0.31 are shown in Fig. 20.As for the densities described in the main article, increasing the effective density leads to an increase in the ellipsoids' local translation order.This can be seen in the first peaks of the pair distribution function g(r) (Fig. 20).However, one also notices that the long range order in the plot of g(r) for the state Γ = 0.31 is weaker than for the nematic states Γ = 0.27, 0.29.This is expected as the ellipsoid alignment in the nematic states enhances the translation order while in the liquid glass there is no long range alignment.Plotting the structure factor S(q) (Fig. 20), we find that the small bump at q < 5 is lower in the liquid glass than in the nematic states.Also, the amplitude of the most intense peak in the S(q) plot is different from density to density.For the nematic phase, the first peak increases with the density but for the liquid glass the first peak is smaller than the first peak in the other nematic states.The fact that the nematic alignments do not extend for long distances in the liquid glass (see Fig. 23) makes the translation order smaller than the one in the nematic states at small wave vectors.
Functions giving insight in the orientation structures for the three states are shown in Fig. 21.The orientation pair distribution function G 2 (r) indicates the nematic order for the states Γ = 0.27, 0.29, 0.31.In the liquid-glass state Γ = 0.31, the drop of the function as the center to center distance r increases is apparent.For the two nematic states Γ = 0.27, 0.29 the nematic order persists and the peaks in the plots of these two states show strong nematic order.For large r, G 2 (r) corresponds to the quadratic value of the scalar nematic order parameter S. Here, we take the nematic order parameter that corresponds to a cluster of size 34 ellipsoids (see Fig. 23) since equilibrating the whole nematic states was too time consuming.For the nematic states, we notice that the peaks of G 2 (r) for the state Γ = 0.29 are higher than the peaks for the state Γ = 0.27 at all length scales.This indicates that the nematic order becomes more enhanced as the nematic density increases.The second orientation structure function is the orientation structure factor S 20 (q) given in Fig. 21.Differences between the three states in the S 20 (q) plot appear at small wave vector values around q < 4. One notices that when increasing the density, the values of S 20 (q) at very small wave vectors become larger.At wave vectors q < 1, the S 20 (q) values are highest for the liquid glass state at Γ = 0.31.This indicates that the fluctuations of the local orientation microscopic densities extend to the whole simulation box and a phase transition takes place in which the global nematic order is destroyed.
The dynamics of the three states Γ = 0.27, 0.29, 0.31 can be inferred from the translation correlation functions F s (q = 2.7, t), S(q = 2.7, t), S 20 (q = 2.7, t) and from the orientation correlation function L 2 (t) in Fig. 22.The wave vector q = 2.7 is close to the first peak in the plots of S 20 (q) in Fig. 21.The three translation correlation functions relax to equilibrium.They decay at longer times when the density increases.For the orientation correlation function L 2 (t), the nematic states decay to the quadratic values of their scalar nematic order parameters S corresponding to the values of S computed for a cluster size of 34 ellipsoids shown in Fig. 23.The reason for not choosing the cluster with the largest number of ellipsoids is that equilibrating the whole nematic samples at such high densities takes a lot of time.On the other hand, the function L 2 (t) for the state Γ = 0.31 does not show any relaxation to equilibrium even though its scalar nematic order parameter S vanishes for the cluster with the largest number of ellipsoids as shown in Fig. 23.In the nematic states, for all cluster sizes the scalar nematic order parameters S are significant.Therefore, we identify Γ = 0.31 as the first liquid glass state in our Brownian dynamics simulations of soft ellipsoid fluid with aspect FIG.22: The self intermediate scattering function Fs(q = 2.7, t), the dynamic structure factor S(q = 2.7, t), the orientation dynamics structure factor S 20 (q = 2.7, t), and the second order orientation correlation functions L 2 (t) for the soft ellipsoid system for the nematic states Γ = 0.27 and Γ = 0.29 and for the liquid glass state Γ = 0.31.S 2 shown in the plot of L 2 (t) is the quadratic value of the scalar nematic order parameter S when the cluster size is 34 particles (see Fig. 23).ratio η = 3.5.We conclude that the liquid glass transition occurred in the density range 0.29 < Γ < 0.31.Fig. 24 shows the mean squared displacements (MSD) for the isotropic states Γ = 0.16, 0.18, the nematic states Γ = 0.21 − 0.29, and the liquid glass states (Γ = 0.31, 0.33).In all states, the Brownian ellipsoids diffuse identically at very short time scales before collisions start to take place at t ≈ 2 × 10 −4 .At longer times, arrest is seen as a weak plateau for the isotropic states and the nematic states Γ = 0.16 − 0.25, and as a clear plateau in the nematic state Γ = 0.29 as well as in the liquid glass states Γ = 0.31, 0.33.At long times, particle motions again become diffusive for all states.Astonishingly, for large time scales t > 3 × 10 2 the liquid glass states Γ = 0.31, 0.33 show identical diffusion independent of the density.The ellipsoids' centers of masses at different densities in the liquid glass for large timescales move similarly.This reminds of plots of the translation correla- tion functions in the liquid glass in the main article which show that the final decays of the translation correlation functions F s (q, t), S(q, t), S 20 (q, t) at large timescales occur identically and independent of the liquid glass density (cf.Fig. 14 and Fig. 15 in the main article).

FIG. 5 :
FIG. 5: Dynamics of the isotropic states in the simulation and experiment.For Γ = 0.16, τ B = 140 sec.For Γ = 0.18, τ B = 220 sec.The experiment (hard ellipsoid) data are the squares and the circles.The solid and the dashed lines are the simulation (soft ellipsoid) curves.

FIG. 6 :
FIG. 6: Dynamics of the isotropic, nematic, and liquid glass states in the simulation.See the text for more details.The isotropic states are the solid lines.The nematic states are the dashed lines with S 2 added from Fig. 16 (at particle number ≈ 275).The liquid glass state Γ = 0.33 is the solid line marked with the filled circles.

FIG. 8 :
FIG.8: Simulation: The scalar nematic order parameter S(t) for different waiting times for Γ = 0.33 (upper panel) and for different densities for the waiting time t 3 w (lower panel).It is clear that the liquid glass does not show any significant nematic order.

17
FIG. 12:The structural functions, g(r) and S(q), at different liquid glass densities.λ is a scaling parameter which maps the simulation (soft-ellipsoid) curves to the experimental (hardellipsoid) curves.For the hard ellipsoids λ = 1.0 and for the soft ellipsoids λ = 1.04.The vertical dashed lines in S(q) mark the wave vectors used in Fig.15, the vertical solid lines mark the q whose dynamics is analysed according to MCT in Fig. 17.
FIG. 12:The structural functions, g(r) and S(q), at different liquid glass densities.λ is a scaling parameter which maps the simulation (soft-ellipsoid) curves to the experimental (hardellipsoid) curves.For the hard ellipsoids λ = 1.0 and for the soft ellipsoids λ = 1.04.The vertical dashed lines in S(q) mark the wave vectors used in Fig.15, the vertical solid lines mark the q whose dynamics is analysed according to MCT in Fig. 17.

FIG. 18 :
FIG.18:The MCT fit to the orientation correlation functions in the simulation at the density Γ = 0.33 for the order n = 10 and n = 20.

FIG. 19 :
FIG. 19: The MCT translation and rotation non-ergodicity parameters (f c l ) and critical amplitudes (h l ) obtained from the MCT fits to the translation and rotation correlation functions at the density Γ = 0.33 for different wave vectors and different Legendre-polynomial orders, respectively; for the fits see Figs.17and 18.The left panel gives the translation amplitudes and the right panel the rotation amplitudes.

17
FIG. 19: The MCT translation and rotation non-ergodicity parameters (f c l ) and critical amplitudes (h l ) obtained from the MCT fits to the translation and rotation correlation functions at the density Γ = 0.33 for different wave vectors and different Legendre-polynomial orders, respectively; for the fits see Figs.17and 18.The left panel gives the translation amplitudes and the right panel the rotation amplitudes.
FIG. 19: The MCT translation and rotation non-ergodicity parameters (f c l ) and critical amplitudes (h l ) obtained from the MCT fits to the translation and rotation correlation functions at the density Γ = 0.33 for different wave vectors and different Legendre-polynomial orders, respectively; for the fits see Figs.17and 18.The left panel gives the translation amplitudes and the right panel the rotation amplitudes.

FIG. 20 :
FIG. 20:The pair correlation function g(r) and the structure factor S(q) for the soft ellipsoid system for the nematic states Γ = 0.27 and Γ = 0.29 and for the liquid glass state Γ = 0.31.No hard ellipsoid (experiment) data are shown here.

FIG. 21 :
FIG. 21:The orientation pair correlation function G 2 (r) and the orientation structure factor S 20 (q) for the soft ellipsoid system for the nematic states Γ = 0.27 and Γ = 0.29 and for the liquid glass state Γ = 0.31.S 2 shown in the plot of G 2 (r) is the quadratic value of the scalar nematic order parameter when the cluster size is 34 particles (see Fig.23).