Structural Correlations and Lack of Order.

To gain insight into the systems’ behavior, we first extracted the pair correlation functions $$ g\left(r\right)$$ and equilibrium structure factors $$ S\left(q\right)$$ from the 3D particle positions at different $$ \varphi $$ (Fig. 1 C and D). With increasing $$ \varphi $$, we observe the emergence of next-neighbor and second next-neighbor peaks in $$ g\left(r\right)$$ while long-range correlations are absent. At the same time, in $$ S\left(q\right)$$, the inverse peak position $$ 2\pi /q_{max}$$, indicating the average particle separation, becomes smaller with increasing density. All observations are typical for liquid-like structural correlations and show that the systems remained in a disordered state without translational order for all investigated $$ \varphi $$. The excellent coincidence between the measured $$ g\left(r\right)$$ and $$ S\left(q\right)$$ with those obtained by the ED-BD simulations confirms that, in the suspensions, our elliptical colloids interacted like hard particles. A further support for the hard particle potentials of our system comes from the lack of formation of nematic or smectic phases known from charged particles (compare below) (39, 40).

The equilibrium phase diagram of ellipsoids predicts the transition to a nematic state with increasing $$ \varphi $$ (19); thus we determined the orientational order of the whole sample. It is probed by calculating the nematic order parameter $$ S$$ (22), which is the eigenvalue possessing the largest absolute value of the nematic order tensor $$ Q_{ij}$$. The nematic order tensor is given by

Fig. 2.

Nematic order parameter $$ S$$ obtained from the nematic tensor $$ {\mathbf{Q}}_{\mathbf{i}\mathbf{j}}$$. (A) $$ S$$ for different $$ \varphi $$, indicating the absence of a global particle alignment in the experimental system. (B) Time evolution of the nematic order parameter for different $$ \varphi $$, again showing no signs of an evolving directional order.

Dynamical Signatures of Glass Formation.

The absence of structural order even in very dense suspensions hints at the possible formation of a glass. To verify this, it is necessary to investigate the structural relaxation dynamics of the suspensions. The translational and rotational dynamics of the system are contained in the temporal evolution of the self part of the density correlation function

Fig. 3 A, C, and E describes the translation behavior and depicts the experimentally determined $$ F_s\left(q,t\right)$$ values for each $$ \varphi $$. In Fig. 3A, the wavevector is set to $$ q=2.6\hspace{0.17em}\mathrm{\mu }{\mathrm{m}}^{-1}$$ (namely, $$ 2bq=6.42$$) which corresponds to the maximum in the static structure factor and thus to the average particle distance. In Fig. 3C, a larger wavevector, $$ q=3.2\hspace{0.17em}\mathrm{\mu }{\mathrm{m}}^{-1}$$ ($$ 2bq=7.90$$), is selected, where the relaxation is faster and the final decay can be resolved well in fluid states. Correlation functions obtained for different values of $$ q$$ can be superimposed by scaling with $$ \log \left(t{q}^2\right)$$ for low $$ \varphi $$ as is expected for the dynamics of isolated ellipsoids (SI Appendix, Fig. S1). For a liquid state, one expects the translation correlators $$ F_s\left(q,t\right)$$ to decay completely over time. During the measurement time, this is observed for volume fractions up to $$ \varphi =0.46$$. The curves obtained for $$ \varphi =0.55$$ are already stretched and thus show a slowing down of the decay dynamics, but one can still assume their complete decay albeit at times longer than those probed. The states at $$ \varphi \le 0.55$$ thus are fluid-like. By contrast, at the highest volume fraction $$ \varphi =0.57$$, the translational correlator does not decay any more, which is seen by the clear plateau for translational motion typical for a glass state. From these observations, we conclude that, for translational dynamics, a glass transition occurred between $$ \varphi =0.55$$ and $$ \varphi =0.57$$.

Fig. 3.

Temporal correlation functions capturing translational and rotational dynamics for different $$ \mathit{\varphi }$$. A and C depict the self part of the density correlation $$ F_s\left(q,t\right)$$ at wavelengths comparable to the average particle separation, namely, $$ 2bq=6.42$$ (A) and $$ 2bq=7.90$$ (C), respectively. B and D show the orientational correlation function $$ L_n\left(t\right)$$ for $$ n=2$$ and $$ n=4$$; the legend in B gives the packing fractions for A–D. Close to the glass transition, the decay of the curves is too slow to be captured within the measurement times. C and D include fits to the correlators in fluid states using the Kohlrausch function Eq. 4. E and F show the relaxation times obtained from the fits. The colored regions mark the glass transitions $$ {\varphi }_c^t$$ and $$ {\varphi }_c^r$$ which were obtained by the MCT glass stability analysis (see Materials and Methods).

Fig. 3 B and D depicts the temporal relaxation behavior of the rotation correlators $$ L_n\left(t\right)$$ for $$ n=2$$ and $$ n=4$$, respectively. A comparison with their translational counterparts, illustrated in Fig. 3 A and C, clearly shows the differences in the translational and rotational relaxation of the suspensions. For volume fractions $$ \varphi \le 0.46$$, a complete relaxation of orientational correlation is observed even if the curves for $$ \varphi =0.46$$ are already stretched significantly. However, for the two highest $$ \varphi $$, a clear plateau indicative of incomplete relaxation appears in both orientational correlation functions. These frozen orientational correlations are clear signs of glass-like behavior. Plateaus in the orientational correlation functions could also arise from nematic ordering. The observed values for the plateau heights, however, are too large to be compatible with the negligible nematic order of our samples. In nematic states, the relation $$ L_2\left(t\to \infty \right)=S^2$$ holds (41), and the small values of the order parameter $$ S$$ from Fig. 2 are incompatible with the high amplitudes of frozen-in orientational correlations. We thus conclude that the data show a glass transition in the rotational dynamics occurring between $$ \varphi =0.46$$ and $$ \varphi =0.55$$.

The inspection of Fig. 3 leads to the conclusion that two different glass transitions, one in the orientational and a second one in the translational motion, exist. To gain more insight into this phenomenon, we turn to a quantitative analysis. In fully relaxing systems, the long-time decay, often termed $$ \alpha $$-relaxation, can be fitted with a Kohlrausch function (42),

A quantitative determination of the $$ \varphi $$ at which the glass transition for translation $$ {\varphi }_c^t$$ and rotation $$ {\varphi }_c^r$$ occurs is obtained by the glass stability analysis of MCT (43) (compare Materials and Methods). Performing such an analysis, we find that satisfactory fits to the data (Fig. 3 E and F) can only be obtained assuming two different glass transitions, $$ {\varphi }_c^t=0.56$$ and $$ {\varphi }_c^r=0.50$$. Therefore, also, this analysis shows the existence of two separate glass transitions for rotation and translation in the experiment.

In summary, the analysis of the systems’ dynamics reveals that we observed a density region where orientational motion in the sample was frozen while translational motion persisted. As this state lacks global nematic order, it is properly described as a liquid glass (21).

ED-BD Simulations Pointing to Nematic Correlations.

The interpretation of experimental results obtained from colloidal suspensions often is complicated by the fact that the particles possess a certain polydispersity and might carry residual charges. Since these problems are nonexistent in simulations, the latter are an important tool even if the number of particles which can be studied in this manner is typically rather small. To understand the origin of the formation of liquid glass in our system, we performed additional ED-BD simulations. We found that nematic order sets in quickly using monodisperse ellipsoids and periodic boundary conditions (PBCs) in accordance with previous simulation studies (22, 44). Therefore, we tested different scenarios that could explain the formation of the experimentally observed glass states. First, since the experimental system is inherently polydisperse, we included polydispersity in our simulations. However, we found that a polydispersity of 5 to 15%, well above that of the experimental particle system used, did not affect the structure and dynamics of the system. Second, in order to prevent the formation of nematic order, we introduced rough walls in our simulations. The approach of confining glass-forming liquids has been used experimentally (45) and in simulations (46) for studies of polydisperse spheres. Indeed, we find that nematic order is suppressed in systems with rough walls. Results for the dynamics of this system are shown in Fig. 4 for $$ F_s\left(q,t\right)$$ and $$ L_4\left(t\right)$$. Here, $$ L_4$$ is shown, as it exhibits glassy dynamics more clearly than $$ L_2$$. The lines represent the MCT fitting curves obtained from Eqs. 6 and 7. Up to $$ \varphi =0.53$$, both correlation functions decayed to fluid states, while, at $$ \varphi =0.55$$ and $$ \varphi =0.57$$, the system was already very close to the glass phase. However, in contrast to the experiments, the simulations still result in a single value for the glass transition of $$ {\varphi }_c^t={\varphi }_c^r\approx 0.57$$.

Fig. 4.

Glass transition analysis of simulation data. (A) Translational $$ F_s\left(q,t\right)$$ for $$ 2qb=6.42$$ and (B) orientational $$ L_4\left(t\right)$$ correlation functions from ED-BD simulations, where long-range nematic fluctuations are suppressed by rough walls. An MCT glass stability analysis (black lines for the two highest densities; compare Materials and Methods) finds a single glass transition density in this case. The value $$ {\varphi }_c\approx 0.57$$ is also supported by a power law fit (with $$ \gamma =2.5$$) to the Kohlrausch relaxation times for both translation and rotation.

While the latter discrepancy between the experimental system and simulations appears dissatisfying, it allows us to formulate a working hypothesis on the prerequisites for the formation of a glass and a liquid glass, respectively: It is known, from experiments on suspensions of spherical colloids, that the glass transition at the higher $$ {\varphi }_c^t$$ results from isotropic caging (42). This is accessible by simulation if nematic order is suppressed by the introduction of rough walls. Under this condition, the build-up of local neighbor shells leading to caging is possible already in small systems (24). The growing fluid structure is reflected in the evolution of the structure factor $$ S\left(q\right)$$ shown in Fig. 1D. By contrast, the formation of a liquid glass above $$ {\varphi }_c^r$$ and below $$ {\varphi }_c^t$$ requires long-range correlations, which cannot be captured in the simulation box sizes currently accessible. This is especially true for shape anisotropic particles, for which it has been shown that simulations are susceptible to finite size effects (47). Solvent-mediated interactions, which were neglected in the simulations, could also play a role, but have been shown to be of minor importance in the structural dynamics of spherical colloids (48). We therefore conjecture that the discrepancy between experiment and simulations reveals that the nature of the observed liquid glass state is relying on long-range correlations.

Nematic Precursor Analysis.

To identify these long-range correlations and to test our hypothesis, we analyzed the spatial correlations in the samples. A close inspection of the imaging data reveals clusters of similarly oriented particles in the system. Two particles are assigned to the same precursor cluster, if they are next neighbors sharing a face of their individual Voronoi cells and have an angular difference in orientation which is less than $$ \mathrm{\Delta }\alpha \le 20°$$. The choice of $$ \mathrm{\Delta }\alpha $$ is based on the half-width of the peak in the probability distribution $$ P\left(\mathrm{\Delta }\theta \right)$$ of the orientation of all ellipsoids in systems at high $$ \varphi $$ (Fig. 5D). Clusters are defined as nematic precursors, if they contain more than 30 particles. The choice of 30 particles is based on the largest cluster sizes found for $$ \varphi =0.40$$, which corresponds to an isotropic sample. Since the clusters do not form nematic regions (Fig. 5H), we term them nematic precursors. Fig. 5 A–C depicts the 20 largest nematic precursors for the three highest $$ \varphi $$ studied. In general, we find that the number of particles in nematic precursors increases with growing $$ \varphi $$, as shown in Fig. 5G. Also, the cluster size distributions become broader with increasing packing fraction. Our findings show that, with $$ \varphi $$ increasing to densities where we find glassy states, more and more particles are found in local structures with a peculiar order. In these, no nematic order is detected even on small length scales, because other particles and clusters intersect the nematic precursors. This is exemplified in Fig. 5E, where the green particles belong to a nematic precursor whereas the blue particles are disordered surrounding particles intersecting the cluster. A local analysis of the nematic order shows this quantitatively. To this end, we divided the sample into equal subboxes and calculated the maximal $$ S_{max}$$ of the order tensor Eq. 1 per subbox. We started with subboxes containing an average of $$ \sim 10$$ particles and gradually increased their size up to $$ \sim 10^3$$ particles. The mean and variance of $$ S_{max}$$ vs. the number of particles contained in the subboxes are plotted in Fig. 5H. The results show a quick drop of the local order parameter $$ S_{max}$$ to low values with increasing particle number. Uniform order does not extend to more than roughly 20 particles. In the vicinity of the glass transition, orientational correlations exist among many particles (Fig. 5G), but still the nematic order parameter remains negligible (Fig. 5H). Nematic order formation is found only for particles with larger aspect ratios (49). The existence of the nematic precursors is also visible in the orientational structure factor, $$ S_{20}\left(q\right)$$, which reflects the density fluctuations of quadrupolar symmetry, $$ {\rho }_{l=2,m=0}\left(\mathbf{\text{q}}\right)$$, proportional to the spherical harmonic for $$ l=2$$ and $$ m=0$$ (21) (SI Appendix). It is sensitive to the growth of nematic structures, which shows up as a peak for small wavevectors, and has been identified as correlation driving the glass transition at $$ {\varphi }_c^r$$ in MCT (21). Fig. 5F (note the log-linear plotting) shows that the large-$$ q$$ peak in $$ S_{20}\left(q\right)$$, which records the orientational alignment of neighboring particles, remains rather unchanged while, upon increasing $$ \varphi $$, the nematic precursors cause a growing peak at small wavevectors. Its width measures the average size of the clusters, which exceeds the average particle separation by roughly a decade. We have carefully evaluated other possible causes for the appearance of long-wavelength fluctuations and could exclude them for the following reasons: 1) Repulsive or attractive forces between the particles would cause differences in the pair correlation and structure factor data derived from experiment and simulation, which is not the case (Fig. 1). 2) The effect of polydispersity was found to be negligible in additional simulations. 3) Scanning electron microscope images (Fig. 1A) show the particles to be smooth on a scale of tens of nanometers, which rules out surface roughness. 4) Thermal fluctuations were fully included in the simulations, and 5) hydrodynamic interactions are considered insignificant at these volume fractions (48). 6) Mermin–Wagner fluctuations, which were suggested to affect glass formation in 2D ellipsoids (50), do not exist in 3D systems. Our data therefore suggest that nematic precursors are the hallmark of the long-wavelength fluctuations observed in bulk liquid glass.

Fig. 5.

Analysis of nematic precursors in the experimental system at high densities. (A–C) A 3D representation of the 20 largest clusters with similar orientation, defined as described in the text, for different $$ \varphi $$ showing the existence of nematic precursors which are intercepting each other. Box size is $$ 140\times 140\times 40\hspace{0.17em}\mathrm{\mu }$$m each. All particles are scaled by a factor of 0.8 for better visibility. (D) Probability distribution $$ P\left(\theta \right)$$ of the orientation of all ellipsoids in the system for different $$ \varphi $$. For increasing $$ \varphi $$, a favored orientation is emerging, indicating nematic precursors but no apparent nematic order. (E) The green particles belong to one cluster with similar orientation, while the surrounding blue particles are in an amorphous arrangement or belong to clusters of different orientations. The hindrance between the different precursors prevents global nematic order. Box size is $$ 25\times 25\times 25\hspace{0.17em}\mathrm{\mu }$$m. (F) Orientational structure factor $$ S_{20}\left(q\right)$$ obtained from spatially correlating the orientations of all particles. Its dominating small-$$ q$$ peak records the growth of the nematic precursors. Inset shows zoom in on low $$ q$$ values. (G) Normalized cluster size $$ N_c/N_{tot}$$ for high $$ \varphi $$, where $$ N_c$$ is the number of particles within a cluster and $$ N_{tot}$$ the total number of particles within the image volume. The cluster size is increasing for higher $$ \varphi $$. (H) Dividing the sample into equal subboxes, the mean and variance of 1) nematic order parameter $$ S_{max}$$ and 2) number of particles per subbox are shown. The 3D liquid glass remains isotropic even for small subboxes. In 2D glasses of ellipsoidal colloids (31), nematic domains containing close to 100 particles exist. The dashed-dotted line marks the isotropic–nematic transition; note the different axes for 3D and 2D.

The finding of a liquid glass state and of nematic precursors in a system of colloidal ellipsoids in three dimensions is in stark contrast to experimental and theoretical results for colloidal ellipsoids confined in two dimensions (31, 51). While, for an aspect ratio of $$ a/b=3.5$$, the 2D studies also show the appearance of a state in which rotational dynamics is frozen and translational mobility persists, the structural order of this state is very different from the 3D case, as, in two dimensions, domains with high nematic order are found. Repeating the local analysis of the nematic order parameter with 2D data from ref. 31, Fig. 5H shows that homogeneously aligned nematic subboxes can comprise up to a hundred particles in 2D glasses. For the comparison in Fig. 5H, we used $$ S_{max}=0.3$$ and 0.5 as a criterion for the isotropic–nematic transition in three dimensions and two dimensions, respectively (22, 52). The analysis reveals that pseudonematic domains exist in two dimensions; the local structure in three dimensions, however, remains isotropic even at the lengthscale of the smallest subbox.

Experimental Methods.

In a first step, spherical PMMA/PMMA core–shell particles dyed with different fluorophores in the cores and the shells were synthesized according to a route described in refs. 35 and 55. As fluorophores for the cores, we used the Bodipy dye ((4,4-difluoro-8-(4-methacrylatephenyl)-3,5-bis-(4-methoxyphenyl)-4-bora-3a,4a-diaza-s-indacene), which was synthesized according to ref. 56. The dye Quasar 670 (LGC Biosearch Technologies) served as a marker for the shells and was used in the form of a 4-vinylbenzylester. For core particle synthesis, the dye was linked to methacrylic acid (Fluka) and copolymerized. The cores were cross-linked using ethylene glycol dimethylacrylate. The resulting particle cores had a diameter of 1$$ \hspace{0.17em}\mathrm{\mu }$$m. Several seeded dispersion polymerization steps were then used to grow a PMMA shell of the desired thickness onto these cores. The PMMA in the shell material was not cross-linked.

In a second step, the particles were thermomechanically stretched as reported by Keville et al. (17), with modifications for the PMMA/PMMA core–shell particles as described in refs. 34, 37, and 57. After the stretching, the ellipsoidal particles were restabilized with polyhydroxystearic acid and dispersed in a mixture of 85:15 (w/w) cyclohexylbromide (CHB) and cis-decalin to match both the density and refractive index of the PMMA colloids (36). Since PMMA particles are known to charge in organic solvents, tetrabutylamoniumbromide was added to screen the charges.

The resulting suspensions were placed into a home-built sample chamber. Samples were studied using a confocal laser scanning microscope (TCSP5, Leica Microsystems) equipped with a resonant scanner (8 kHz, bidirectional mode) and a glycerol objective (63×, NA 1.30). An argon laser (wavelength $$ \lambda =514.5$$ nm) and a helium neon laser ($$ \lambda =632.8$$ nm) were used for excitation of the dyes. The temperature of the whole microscope system was stabilized at $$ 23.0\pm 0.1\hspace{0.17em}°$$C (Ludin Cube 2, Life Imaging Services), since, at this temperature, the density match between particles and solvent was optimal. Volume fractions up to $$ \varphi =0.46$$ were prepared by centrifugation at higher temperature and the removal of excess solvent. The highest $$ \varphi $$ were obtained by slightly raising the temperature in a dense sample directly in the measurement chamber for a few days to achieve density mismatch. This led to the slow enrichment of particles at the bottom of the sample.

The temperature was then changed back to $$ 23\hspace{0.17em}°$$C. After placing a sample onto the microscope, it was allowed to equilibrate for 15 h to minimize drift effects. Remaining drift, which, in all cases, was small compared to the particle movement, was corrected with a routine applied after detection and tracking. To minimize wall effects, all measurements were performed in a depth of 30 $$ \mathrm{\mu }$$m to 100 $$ \mathrm{\mu }$$m into the sample. Positions and orientations of all particles in one image volume (144 $$ \mathrm{\mu }$$m $$ \times 144\hspace{0.17em}\mathrm{\mu }$$m $$ \times 40\hspace{0.17em}\mathrm{\mu }$$m for structural analysis and 144 $$ \mathrm{\mu }$$m $$ \times 36\hspace{0.17em}\mathrm{\mu }$$m $$ \times 20\hspace{0.17em}\mathrm{\mu }$$m for measurements on the dynamics) containing more than 6,000 particles were identified and linked to trajectories using self-written and tested Matlab algorithms. The core–shell structure facilitates tracking of the ellipsoid position and orientation with accuracies of 60 nm and $$ 5°$$, respectively (35) (Fig. 1B). Experimental timescales were limited by the bleaching of the fluorophores after long laser exposure.

Volume fractions are the pivotal variable influencing the structural dynamics of the colloidal suspensions. For their determination, particles were detected and counted in each measured sample, and the related Voronoi volume was calculated to obtain the exact number fraction $$ \rho $$. Multiplying this with the volume of a particle gives $$ \varphi $$,

Glass Stability Analysis.

The values for the glass transitions for translation $$ {\varphi }_c^t$$ and rotation $$ {\varphi }_c^r$$ are found by the glass stability analysis of MCT in which the collection of all measured correlation functions is used. Close to the glass transition $$ {\varphi }_c$$, a critical decay onto the plateau occurs in a correlation function $$ {\mathrm{\Phi }}_j\left(t\right)$$ (e.g., the translation or the rotation correlator functions $$ F_s\left(q,t\right)$$ and $$ L_n\left(t\right)$$, respectively). Using first-order correction terms (43), the decay is given by

Below $$ {\varphi }_c$$, the late plateau and the decay from the plateau is given by the von Schweidler law,

Fig. 6.

Glass stability analysis of (A) translational and (B) rotational relaxations: Data from experiment (colored lines) and fits (solid lines) according to Eqs. 6 and 7. Two different glass transitions for translation $$ {\varphi }_c^t=0.56$$ and rotation $$ {\varphi }_c^r=0.50$$ are necessary to fit the data. In C, the gray dash-dotted lines indicate the result, when an identical transition density is assumed for translation and rotation. For this condition, poor agreement with the data at $$ \varphi =0.46$$ is observed.

Simulation Details.

We performed ED-BD simulations (20, 38) that model overdamped Brownian dynamics of hard ellipsoids. The system consisted of $$ N$$ nonoverlapping ellipsoids in a simulation box of length $$ L$$ and volume $$ V=L^3$$ that was varied depending on the desired $$ \varphi $$. Similar to De Michele et al. (22), we set the moment of inertia to $$ I=2m{b}^2/5$$. For volume fractions of $$ \varphi =0.05,0.1,0.2$$, ellipsoids were placed randomly in the simulation box, while, for volume fractions between $$ \varphi =0.4$$ and $$ \varphi =0.57$$, the Lubachevsky–Stillinger (58) technique was used. All systems were equilibrated by Newtonian dynamics, allowing translational and rotational correlation functions to decay to zero. Two types of boundary conditions were modeled. The first case corresponds to systems with PBCs with $$ N=504$$. The second case corresponds to systems confined within rough walls. For the confined systems, we started from a random configuration of 1,231 ellipsoids and pinned the outermost particles in the range $$ 0<d<\mathrm{\Delta }$$ and $$ L-\mathrm{\Delta }<d<L$$ to serve as the wall particles. The remaining 500 particles were free to move. The ED-BD model was first benchmarked for a dilute system of $$ \varphi =0.05$$, where we reproduced the ratio of translational and rotational diffusion coefficients found in the experiments: $$ D_{\text{rot}}{b}^2/D_{\text{trans}}=1.25$$, setting the ratio of Brownian time steps $$ \mathrm{\Delta }{t}_{\text{rot}}/\mathrm{\Delta }{t}_{\text{trans}}=0.331$$ (with $$ \mathrm{\Delta }{t}_{\text{rot}}=2I/k_{B}T{D}_{\text{rot}}=0.0265$$ short enough to give angular diffusion). This ratio was then applied to all $$ \varphi $$. At all densities, the timescale, when matching to the experiments, is taken from $$ D_{\text{trans}}\left(\varphi \right)$$. All presented results were averaged from 20 initial configurations of the ellipsoids.