You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'We numerically simulate the uniform athermal shearing of bidisperse, frictionless, two dimensional spherocylinders and three dimensional prolate ellipsoids. We focus on the orientational ordering of particles as an asphericity parameter $\alpha\to 0$ and particles approach spherical. We find that the nematic order parameter $S_2$ is non-monotonic in the packing fraction $\phi$, and that as $\alpha\to 0$ $S_{2}$ stays finite at jamming and above. The approach to spherical particles thus appears to be singular. We also find that sheared particles continue to rotate above jamming, and that particle contacts preferentially lie along the narrowest width of the particles, even as $\alpha\to 0$.'
author:
- Theodore Marschall
- 'Yann-Edwin Keta'
- Peter Olsson
- 'S. Teitel'
title: 'Orientational Ordering in Athermally Sheared, Aspherical, Frictionless Particles'
---

Models of athermal ($T=0$) granular materials have often focused on the simplest case of spherical particles. Recently, however, more attention has been paid to the case of elongated particles with lower rotational symmetry [@Borzsonyi.Soft.2013]. The question then arises whether such elongated particles will orientationally order as the particle density increases, so as to pack more efficiently. While elongated particles in thermal equilibrium are known to have a nematic orientational ordering transition [@Onsager; @Bolhuis], recent works have found that such particles do [*not*]{} orientationally order upon athermal isotropic compression [@Donev.Science.2004; @Man.PRL.2005; @Sacanna.JPhysC.2007; @Zhao; @Marschall].

Orientational ordering is, however, found when elongated particles are placed in an athermal uniform shear flow. In this case, drag forces between the particle and the flow will cause the particle to tumble [@Jeffery]. If the particle is asymmetrical, with unequal eigenvalues of its moment of inertia tensor, tumbling will have a non-uniform rotational velocity; the torque from drag forces will vary with the particle’s orientation, and the particle will on average align with the flow direction. For a finite density of colliding particles, nematic ordering remains but the nematic director becomes oriented at a finite angle with respect to the flow direction [@Campbell; @Guo1; @Guo2; @Borzsonyi1; @Borzsonyi2; @Wegner; @Wegner2; @Nagy; @Trulsson].

Here we investigate the nematic ordering of frictionless, aspherically shaped particles, athermally sheared at constant strain rate $\dot\gamma$, putting the system into a steady state of simple shear flow. We consider behavior as an asphericity parameter $\alpha\to 0$, and the particles approach spherical. We find the surprising result that a [*finite*]{} nematic ordering persists even as $\alpha\to 0$, suggesting that the shear driven jamming of aspherical particles has a singular limit as $\alpha\to 0$. Since most particles in nature are not truly spherical, our result may have broad implications for granular shear flows.

[**Models:**]{} We consider two different numerical models: (i) spherocylinders in two dimensions (2D), and (ii) prolate ellipsoids in three dimensions (3D). In both cases we take a bidisperse distribution of particle sizes, with equal numbers of big and small particles. We use soft-core particles with a one-sided harmonic elastic repulsion. The system length is $\mathcal{L}$ in all directions, with periodic boundary conditions along the flow direction $\mathbf{\hat x}$, and Lees-Edwards boundary conditions [@LeesEdwards] with a uniform strain rate $\dot\gamma$ in the transverse direction $\mathbf{\hat y}$. In 3D we take periodic boundary conditions along $\mathbf{\hat z}$. The particle packing fraction is $\phi=\sum_i v_i/\mathcal{V}$, with $v_i$ the volume of particle $i$ and $\mathcal{V}=\mathcal{L}^d$ the system volume ($d=2$ or $3$ for 2D and 3D respectively).

[**2D Spherocylinders:**]{} A 2D spherocylinder consists of a rectangle of length $L$, with two semi-circular end caps of diameter $D$ (see inset to Fig. \[f5\]a). We define the asphericity parameter $\alpha=L/D$. Big and small particles have equal $\alpha$, with $D_b/D_s=1.4$. Taking the “spine" of the spherocylinder as the line bisecting the rectangle parallel to its length $L$, we define $r_{ij}$ as the shortest distance between the spines of spherocylinders $i$ and $j$ and $d_{ij}=(D_i+D_j)/2$. Two spherocylinders are in contact whenever $r_{ij}<d_{ij}$, in which case the elastic interaction is $U^\mathrm{el}=(k_e/2)(1-r_{ij}/d_{ij})^2$ and the repulsive elastic force on $i$ is $\mathbf{F}^\mathrm{el}_{ij}=(k_e/d_{ij})(1-r_{ij}/d_{ij})\mathbf{\hat n}_{ij}$, with $\mathbf{\hat n}_{ij}$ the unit vector pointing normally inwards to particle $i$ at the point of contact with $j$ [@Marschall; @SupMat].

Our dynamics is the mean-field Durian model for foams [@Durian], generalized to non-spherical particles. The dissipative force on a spherocylinder is a Stokes drag between the particle and a uniform background shear flow, $\mathbf{F}_i^\mathrm{dis}=-k_d v_i(\dot{\mathbf{r}}_i-y_i\dot\gamma\mathbf{\hat x})$, with $\mathbf{r}_i=(x_i,y_i)$ the center of mass of spherocylinder $i$, $\dot{\mathbf{r}}_i$ the center of mass velocity, and $k_d$ the viscous coupling. We use overdamped dynamics $\mathbf{F}_i^\mathrm{dis}+\sum_j\mathbf{F}_{ij}^\mathrm{el}=0$, where the sum is over all particles $j$ in contact with $i$.

The elastic and dissipative forces produce torques on the spherocylinders. The elastic torque on particle $i$ due to contact with $j$ is, $\boldsymbol{\tau}_{ij}^\mathrm{el}=\mathbf{\hat z}\tau_{ij}^\mathrm{el}=\mathbf{s}_{ij}\times \mathbf{F}_{ij}^\mathrm{el}$, where $\mathbf{s}_{ij}$ is the moment arm from the center of mass of $i$ to its point of contact with $j$. A dissipative torque arises from the variation of the background shear flow velocity over the spatial extent of the particle [@Marschall2]. Integrating over particle area gives $\tau_i^\mathrm{dis}=-k_dv_i I_i[\dot\theta_i+\dot\gamma f(\theta_i)]$, where $\theta_i$ is the angle of the spine with respect to the flow direction $\mathbf{\hat x}$, and $f(\theta)=[1-C\cos2\theta]/2$. The overdamped $\tau_i^\mathrm{dis}+\sum_j \tau_{ij}^\mathrm{el}=0$ determines the particle rotation. Here $I_i$ is the sum of the two eigenvalues of the moment of inertia tensor, and $C=\Delta I_i/I_i$, with $\Delta I_i$ the difference between the two eigenvalues. For spherocylinders, $I_i=(D_i/2)^2(3\pi+24\alpha+6\pi\alpha^2+8\alpha^3)/(6\pi+24\alpha)$. For a circle, $\Delta I=0$, and so in the absence of collisions $\dot\theta/\dot\gamma=-1/2$. We take as unit of length $D_s=1$, unit of energy $k_e=1$, and unit of time $t_0=D_s^2k_d/k_e=1$. We integrate the equations of motion using the Heun method with step size $\Delta t/t_0=0.02$. We use $N=1024$ particles.

[**3D Prolate Ellipsoids:**]{} We take prolate ellipsoids of revolution with major axis length $a_1$ and minor axes length $a_2$. The asphericity is $\alpha=a_1/a_2-1$. Big and small particles have equal $\alpha$, with $a_{1b}/a_{1s}=1.4$. When two ellipsoids $i$ and $j$ overlap, we define a scale factor $\delta_{ij}<1$ such that the particles just barely touch when their axes are rescaled by $\delta_{ij}$, keeping the center of mass positions fixed [@SupMat]. The elastic interaction is then $U^\mathrm{el}=(k_e/2)(1-\delta_{ij})^2$, and the repulsive elastic force on $i$ is $\mathbf{F}_{ij}^\mathrm{el}=k_e \delta_{ij} (1-\delta_{ij})\mathbf{\hat n}_{ij}/[(\mathbf{r}_i-\mathbf{r}_j)\cdot\mathbf{\hat n}_{ij}]$, with $\mathbf{r}_i$ the center of mass of ellipsoid $i$ and $\mathbf{\hat n}_{ij}$ the unit vector pointing normally inwards to particle $i$ at the point of contact with $j$.

We take a purely collisional dynamics. The dissipative force on ellipsoid $i$ is due to contact with $j$ and is proportional to the difference in particle velocities at their point of contact, $\mathbf{F}_{ij}^\mathrm{dis}=-k_d (\dot{\mathbf{r}}_i+\boldsymbol{\omega}_i\times \mathbf{s}_{ij}-\dot{\mathbf{r}}_j-\boldsymbol{\omega}_j\times\mathbf{s}_{ji})$, with $\dot{\mathbf{r}}_i$ the center of mass velocity, $\boldsymbol{\omega}_i$ the angular velocity about the center of mass, and $\mathbf{s}_{ij}$ the moment arm from the center of $i$ to the point of contact with $j$ [@CDrot]. We use Newton’s equation of motion, $m_i\ddot{\mathbf{r}}_i=\sum_j[\mathbf{F}_{ij}^\mathrm{dis}+\mathbf{F}_{ij}^\mathrm{el}]$, where the sum is over all particles $j$ in contact with $i$, and the mass $m_i$ is taken proportional to the particle volume $v_i$. The rotation of particle $i$ is governed by, $\mathbf{I}_i\cdot\dot{\boldsymbol{\omega}}_i=\sum_j\mathbf{s}_{ij}\times [\mathbf{F}_{ij}^\mathrm{dis}+\mathbf{F}_{ij}^\mathrm{el}$\], where $\mathbf{I}_i$ is the moment of inertia tensor of $i$.

We take as unit of length $D_s\equiv\sqrt[3]{a_{1s}a_{2s}^2}=1$, unit of energy $k_e=1$, unit of mass $m_s=1$ and unit of time $t_0=D_s\sqrt{m_s/k_e}=1$. Collision elasticity is measured by $Q=\sqrt{m_s k_e}/(k_d D_s)=2$, which would be the quality factor of a corresponding damped oscillator. We integrate the equations of motion using a modified velocity Verlet algorithm [@CDrot] with step size $\Delta t/t_0=0.05$. We use $N=1024$ particles.

[**Results:**]{} In this work we focus on the orientational order and tumbling of particles, rather than rheology. To measure nematic ordering we compute the tensor [@Borzsonyi1], $$\langle T_{\mu\nu}\rangle=\left\langle\frac{d}{(d-1)N}\sum_{i=1}^N\left[\mathbf{\hat\ell}_{i\mu}\mathbf{\hat\ell}_{i\nu}-\frac{1}{d}\delta_{\mu\nu}\right]\right\rangle,
\label{eT}$$ where $\boldsymbol{\hat\ell}_i$ is a unit vector along the spine of the spherocylinder or the major axis of the ellipsoid, $\mu$ and $\nu$ denote spatial components, $d=2,3$ is the spatial dimension, and $\langle\dots\rangle$ denotes an average over configurations in the sheared ensemble. The largest eigenvalue of $\langle T_{\mu\nu}\rangle$ is the magnitude of the nematic order parameter $S_2$. The corresponding eigenvector $
\boldsymbol{\hat\ell}_2$ gives the orientation of the nematic director, which by symmetry lies in the $xy$ plane; $\theta_2$ is the angle of $
\boldsymbol{\hat\ell}_2$ with respect to the flow direction $\mathbf{\hat x}$, and $\mathbf{S}_2=S_2 \boldsymbol{\hat\ell}_2$.

![Nematic order parameter $S_2$ vs packing $\phi$ at different shear strain rates $\dot\gamma$. (a) 2D spherocylinders at asphericity $\alpha=0.01$, (b) 3D ellipsoids at $\alpha=0.02$. Vertical dashed lines locate the jamming transition of $\alpha=0$ spherical particles, $\phi_J^{(0)}=0.8433$ for 2D [@OT1; @OT2; @OT3] and 0.649 for 3D [@Olsson]. []{data-label="f1"}](f1){width="3.5in"}

In Fig. \[f1\] we plot $S_2$ vs $\phi$ for particles of fixed small asphericity $\alpha$, at different shear rates $\dot\gamma$: (a) 2D spherocylinders at $\alpha=0.01$, (b) 3D ellipsoids at $\alpha=0.02$. Both cases show similar behavior. In contrast to previous works [@Campbell; @Guo1; @Guo2] that found increasing orientational order with increasing $\phi$, here we find a non-monotonic $S_2$ [@Trulsson] with peak value $S_{2\,\mathrm{max}}$ at a $\phi_\mathrm{max}$ slightly below the $\alpha=0$ jamming transition at $\phi_J^{(0)}$. As $\dot\gamma$ decreases, the values of $S_2$ approach a common limiting curve [@Borzsonyi1; @Borzsonyi2]; above $\phi_J^{(0)}$ nematic order $S_2$ stays finite, but there is a stronger $\dot\gamma$ dependence.

In Fig. \[f2\] we plot $S_2$ vs $\phi$ for a range of $\alpha$, showing results for both a smaller $\dot\gamma_1$ (solid symbols) and a larger $\dot\gamma_2$ (open symbols); see Table \[tab1\] for values. In each case $\dot\gamma_1$ and $\dot\gamma_2$ are sufficiently small that $S_2$ shows no noticeable $\dot\gamma$ dependence for $\phi$ up to and slightly beyond the peak at $\phi_\mathrm{max}$, however some small $\dot\gamma$ dependence remains at the highest $\phi$. What is remarkable is that the orientational ordering $S_{2\,\mathrm{max}}$ remains quite sizable even for particles close to spherical with $\alpha=0.001$.

    2D: $\alpha$      $\dot\gamma_1$      $\dot\gamma_2$       3D: $\alpha$      $\dot\gamma_1$      $\dot\gamma_2$
  ----------------- ------------------- ------------------- ------------------ ------------------- -------------------
        0.001        $1\times 10^{-7}$   $4\times10^{-7}$    $\alpha\le 0.02$   $1\times 10^{-7}$   $2\times 10^{-7}$
        0.01         $4\times 10^{-7}$   $1\times10^{-6}$          0.05         $5\times 10^{-7}$   $1\times 10^{-6}$
   $\alpha\ge0.06$   $1\times 10^{-5}$   $4\times 10^{-5}$        $0.2$         $2\times 10^{-6}$   $5\times 10^{-6}$
                                                                   0.7          $5\times 10^{-6}$   $1\times 10^{-5}$

  : Strain rate values used for data in Figs. \[f2\] and \[f3\]

\[tab1\]

![Nematic order parameter $S_2$ for (a) 2D spherocylinders and (b) 3D ellipsoids vs packing $\phi$ for different asphericities $\alpha$, at two different small strain rates $\dot\gamma_1$ (solid symbols) $< \dot\gamma_2$ (open symbols), see Table \[tab1\] for values. Vertical dashed lines locate the jamming $\phi_J^{(0)}$ of spherical particles. []{data-label="f2"}](f2){width="3.5in"}

Fig. \[f2\] shows $S_2$ averaged over the steady state ensemble. We have also computed the instantaneous $S_2(\gamma)$ and $\theta_2(\gamma)$ as functions of the system strain $\gamma=\dot\gamma t$. We find that near and above the peak at $\phi_\mathrm{max}$, $\mathbf{S}_2(\gamma)$ shows random fluctuations about a well defined average; there is no macroscopically coherent tumbling of particles [@SupMat2]. One can still ask if individual particles tumble incoherently [@Borzsonyi2; @Wegner], or whether they are orientationally locked into small fluctuations about the nematic director $\boldsymbol{\hat\ell}_2$. We find the former to be the case.

![Component of average particle angular velocity in the direction of the system vorticity, scaled by strain rate, $-\langle\omega_{zi}\rangle/\dot\gamma$ for (a) 2D spherocylinders and (b) 3D ellipsoids vs packing $\phi$ for different asphericities $\alpha$, at two different small strain rates $\dot\gamma_1$ (solid symbols) $< \dot\gamma_2$ (open symbols), see Table \[tab1\] for values. Vertical dashed lines locate the jamming $\phi_J^{(0)}$ of spherical particles. []{data-label="f3"}](f3){width="3.5in"}

In Fig. \[f3\] we plot the component of the average particle angular velocity in the direction of the system vorticity, scaled by the strain rate, $-\langle\omega_{zi}\rangle/\dot\gamma$; note, $-\langle\omega_{zi}\rangle>0$ indicates clockwise rotation. For 2D spherocylinders, $\omega_{zi}=\dot\theta_i$. In each case we show results at two different strain rates $\dot\gamma_1<\dot\gamma_2$, as in Fig. \[f2\] (see Table \[tab1\]), and find little dependence on $\dot\gamma$ except for the largest $\phi$. Comparing Figs. \[f2\] and \[f3\] we find that rotation velocity is anti-correlated with orientational order; $-\langle\omega_{zi}\rangle/\dot\gamma$ is non-monotonic in $\phi$ and is smallest when $S_2$ is largest. For small but finite $\alpha$, $-\langle\omega_{zi}\rangle/\dot\gamma$ approaches the spherical limit of 1/2 at small $\phi$, but shows a significant dip below 1/2 at $\phi_\mathrm{max}$. For 2D spherocylinders this dip remains sizable even for very small $\alpha=0.001$. For 3D ellipsoids we cannot get accurate results at similar small values of $\alpha$, but Fig. \[f3\]b shows that the trends appear to be the same. We conclude that particles continue to rotate, with finite $\langle\omega_{zi}\rangle/\dot\gamma$, even above jamming.

Returning to the nematic ordering, in Fig. \[f4\]a we plot $S_{2\,\mathrm{max}}$ vs $\alpha$ for both 2D spherocylinders and 3D ellipsoids. Solid lines are fits to the empirical form $S_{2\,\mathrm{max}}=S_0+c\alpha^\beta$, using the five smallest $\alpha$ points. We find $S_0=0.25$ for 2D spherocylinders and $S_0=0.16$ for 3D ellipsoids. If we exclude the data point at the smallest $\alpha=0.001$, then our data would be reasonably fit (dashed lines in Fig. \[f4\]a) by a pure power law with exponent $\approx 0.14$. However, in [@SupMat3] we give detailed tests confirming that our data point at $\alpha=0.001$ is accurate and so should not be excluded.

In Fig. \[f4\]b we plot $\phi_J^{(0)}-\phi_\mathrm{max}$ vs $\alpha$, where $\phi_J^{(0)}$ is the jamming transition for spherical particles. In both 2D and 3D we find $\phi_J^{(0)}-\phi_\mathrm{max}\to 0$ as $\alpha\to 0$, showing that the peak in $S_2$ approaches the jamming transition as $\alpha\to 0$. For 2D spherocylinders we find a power law dependence, $\phi_J^{(0)}-\phi_\mathrm{max}\sim \alpha^\Delta$ with $\Delta\approx 0.43$, as illustrated by the dashed line in the figure. For 3D ellipsoids, our data do not suggest any clear form for the small $\alpha$ behavior. The observations of Figs. \[f2\] and \[f4\] thus lead us to conclude that, even as $\alpha\to 0$ and particles are approaching the spherical limit, a finite nematic ordering $S_2$ exists at the jamming $\phi_J^{(0)}$ and above.

![For 2D spherocylinders and 3D ellipsoids: (a) $S_{2\,\mathrm{max}}$ vs $\alpha$. Solid lines are fits to $S_0+c\alpha^\beta$, using the five smallest $\alpha$ points. Dropping the point at $\alpha=0.001$, dashed lines show power law fits. (b) $\phi_J^{(0)}-\phi_\mathrm{max}$ vs $\alpha$, with $\phi_J^{(0)}$ the $\alpha=0$ jamming point. Solid lines connect the data points; dashed line for the 2D spherocylinders is a power law fit to the five smallest $\alpha$ points. []{data-label="f4"}](f4){width="3.5in"}

To look for a microscopic signature of this singular $\alpha\to 0$ limit, we measure the location on a particle’s surface of the inter-particle contacts. For 2D spherocylinders we define $(r,\vartheta)$ as the radial distance and polar angle of a point on the surface with respect to the center of the particle and the direction of the spine. We define $\mathcal{P}(\vartheta)$ as the probability density per unit surface length to have a contact at $\vartheta$, with normalization $1=\mathcal{A}^{-1}\int_0^{2\pi} d\vartheta\, \sqrt{r^2 + (dr/d\vartheta)^2}\,
\mathcal{P}(\vartheta)$, with $\mathcal{A}$ the perimeter length [@comment]. For 3D ellipsoids, we define $(r,\vartheta,\varphi)$ as the spherical coordinates with respect to the major axis; $\mathcal{P}(\vartheta,\varphi)$ is the probability density per unit surface area to have a contact at $(\vartheta, \varphi)$, with normalization $1=\mathcal{A}^{-1}\int_0^{2\pi}d\varphi\int_0^\pi d\vartheta \sin\vartheta\, r\sqrt{r^2 + (dr/d\vartheta)^2}\,\mathcal{P}(\vartheta,\varphi)$; $\mathcal{A}$ is the surface area. For simplicity we consider $\mathcal{P}(\vartheta)=(2\pi)^{-1}\int_0^{2\pi}d\varphi \mathcal{P}(\vartheta,\varphi)$. For a uniform probability desnity, such as would be for spherical particles, $\mathcal{P}(\vartheta)=1$ in both 2D and 3D.

![Probability $\mathcal{P}(\vartheta)$ for a particle to have a contact at polar angle $\vartheta$ on its surface, for different asphericities $\alpha$ at fixed $\phi$ near $\phi_J^{(0)}$: (a) 2D spherocylinders at $\phi=0.843$ and (b) 3D ellipsoids at $\phi=0.648$, for sufficiently small $\dot\gamma$ that $\mathcal{P}(\vartheta)$ becomes independent of $\dot\gamma$. In (a) the sharp peaks near $\vartheta=\pi/6$ and $5\pi/6$ are shadow effects from particles in contact at $\vartheta=\pi/2$. ](f5){width="3.5in"}

. \[f5\]

In Fig. \[f5\] we plot $\mathcal{P}(\vartheta)$ vs $\vartheta$ for different asphericities $\alpha$ at fixed $\phi$ near $\phi_J^{(0)}$. For each $\alpha$ we use a $\dot\gamma$ sufficiently small that $\mathcal{P}(\vartheta)$ has approached its $\dot\gamma\to 0$ limiting distribution. Unlike the uniform distribution for spheres, we see a sharp peak at $\vartheta=\pi/2$ whose height steadily increases as $\alpha$ decreases. Thus as particles become increasingly spherical, particle contacts increasingly prefer to form along the narrowest width of the particle rather than uniformly over the particle’s surface [@contacts]. The onset of this effect occurs as $\phi$ increases above the jamming $\phi_J$ [@SupMat4]. We note that similar results for $\mathcal{P}(\vartheta)$ have been reported [@Marschall; @OHern] in static, isotropically jammed packings, but in that case there is no nematic ordering and $S_2=0$. One possible explanation for this difference is that it is the anisotropy of the stress in a sheared system, as manifested by directed force chains, that determines a particular direction and gives rise to a non-vanishing $S_2$. Indeed we find that for small $\alpha$ close to and above jamming (but not well below jamming and not for larger $\alpha$), the orientation $\theta_2$ of the nematic director aligns with the minimum stress axis of the stress tensor, which is at $45^\circ$ with respect to the flow direction $\mathbf{\hat x}$.

To examine the role that stress anisotropy plays, we have carried out preliminary simulations of 2D spherocylinders under a [*pure*]{} shear, compressing our system in the $\mathbf{\hat y}$ direction while expanding in the $\mathbf{\hat x}$ direction, both at constant rate $\dot\gamma/2=5\times 10^{-7}$ so as to keep constant area. While simple shear creates a vorticity in the affine velocity field that drives the continuous rotation of individual particles (as in our Fig. \[f3\]), such vorticity is absent in pure shear; we thus find $\langle\dot\theta_i\rangle=0$, the nematic director aligns with the minimal stress axis, and the magnitude $S_2$ is large at small $\phi$, monotonically [*decreasing*]{} as $\phi$ increases. We find $S_2$ from pure shear and simple shear qualitatively agree [*only*]{} when one is close to or above the jamming $\phi_J$, where behavior is likely dominated by extended force chains that restrict particle alignment. For small $\alpha$, pure shear and simple shear differ most at lower $\phi$: for pure shear particles decay to a fixed orientation giving large $S_2$ and $\theta_2=0$, while for simple shear particles continuously rotate averaging out to a small $S_2$; as $\phi$ increases, elastic collisions increase, the rotation slows and becomes more non-uniform, and $S_2$ increases. The non-monotonic behavior of $S_2$ with a peak at $\phi_\mathrm{max}$ is thus a direct consequence of the rotational drive that is present in simple shear but absent in pure shear. See further details in [@SupMat5].

To conclude, we have considered the athermal uniform shearing of bidisperse, aspherical particles in 2D and 3D. A finite particle asphericity $\alpha$ breaks rotational symmetry, and as in earlier works [@Campbell; @Guo1; @Guo2; @Borzsonyi1; @Borzsonyi2; @Wegner; @Wegner2; @Nagy; @Trulsson] we find a finite nematic ordering $S_2$. However one would naively expect that $S_2\to 0$ as the symmetry breaking parameter $\alpha\to 0$. In contrast, here we show that $S_2$ remains finite at jamming and above even as $\alpha\to 0$. This may be viewed in analogy with an Ising model, where the magnetization $m$ stays finite even as the ordering magnetic field $h\to 0$ for $T<T_c$. However there are two significant differences: (i) In the Ising model with $h\to0$, one has $m\to 0$ as $T\to T_c$ from below, while here as $\alpha\to 0$ we find $S_2$ stays finite as $\phi\to \phi_J^{(0)}$ from above; (ii) ordering in the Ising model arises from a microscopic spin-spin interaction that prefers alignment even when $h=0$, while here the microscopic interaction that prefers alignment of the particle major axes would naively seem to vanish as $\alpha\to 0$ and the particles become spherical (though the behavior of $\mathcal{P}(\vartheta)$ suggests that a local ordering interaction may indeed persist even as $\alpha\to 0$).

It would be interesting to see how robust this effect is to the introduction of additional sources of fluctuation, such as a polydispersity in $\alpha$ [@poly], or the presence of thermal effects. We leave such questions to future research.

Simulations were carried out on resources of the Center for Integrated Research Computing at the University of Rochester and of the Swedish National Infrastructure for Computing (SNIC) at HPC2N. This work was supported in part by National Science Foundation Grant No. CBET-1435861.

[99]{}

For a review, see: T. B[ö]{}rzs[ö]{}nyi and R. Stannarius, “Granular materials composed of shape-anisotropic grains," Soft Matter [**9**]{}, 7401 (2013). L. Onsager, “The effects of shape on the interaction of colloid particles," Annals of the New York Academy of Sciences [**51**]{}, 627 (1949). P. Bolhuis and D. Frenkel, “Tracing the phase boundaries of hard spherocylinders," The Journal of Chemical Physics [**106**]{}, 666 (1997). A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin, “Improving the density of jammed disordered packings using ellipsoids," Science [**303**]{}, 990 (2004). W. Man, A. Donev, F. H. Stillinger, M. T. Sullivan, W. B. Russel, D. Heeger, S. Inati, S. Torquato, and P. M. Chaikin, “Experiments on random packings of ellipsoids," Phys. Rev. Lett. [**94**]{}, 198001 (2005). S. Sacanna, L. Rossi, A. Wouterse, and A. P. Philipse, “Observation of a shape-dependent density maximum in random packings and glasses of colloidal silica ellipsoids," J. Phys.: Condens. Matter [**19**]{}, 376108 (2007). J. Zhao, S. Li, R. Zou, and A. Yu, “Dense random packings of spherocylinders," Soft Matter [**8**]{}, 1003 (2012).

T. Marschall and S. Teitel, “Compression-driven jamming of athermal frictionless spherocylinders in two dimensions," Phys. Rev. E [**97**]{}, 012905 (2018). G. B. Jeffery, “The Motion of ellipsoidal particles immersed in a viscous fluid," Proc. R. Soc. A [**102**]{}, 161 (1922). C. S. Campbell, “Elastic granular flows of ellipsoidal particles," Physics of Fluids [**23**]{}, 013306 (2011). Y. Guo, C. Wassgren, W. Ketterhagen, B. Hancock, B. James and J. Curtis, “A Numerical study of granular shear flows of rod-like particles using the discrete element method," J. Fluid. Mech. [**713**]{}, 1 (2012). Y. Guo, C. Wassgren, B. Hancock, W. Ketterhagen and J. Curtis, “Granular shear flows of flat disks and elongated rods without and with friction," Physics of Fluids [**25**]{}, 063304 (2013). T. Börzsönyi, B. Szabó, G. Törös, S. Wegner, J. Török, E. Somfai, T. Bien and R. Stannarius, “Orientational order and alignment of elongated particles induced by shear," Phys. Rev. Lett. [**108**]{}, 228302 (2012). T. Börzsönyi, B. Szabó, S. Wegner, K. Harth, J. Török, E. Somfai, T. Bien and R. Stannarius, “Shear-induced alignment and dynamics of elongated granular particles," Phys. Rev. E [**86**]{}, 051304 (2012). S. Wegner, T. B[ö]{}rzs[ö]{}nyi, To. Bien, G. Rose, and R. Stannarius, “Alignment and dynamics of elongated cylinders under shear," Soft Matter [**8**]{}, 10950 (2012). S. Wegner, R. Stannarius, A. Boese, G. Rose, B. Szab[ó]{}, E. Somfai, and T. B[ö]{}rzs[ö]{}nyi, “Effects of grain shape on packing and dilatancy of sheared granular materials," Soft Matter [**10**]{}, 5157 (2014). A similar non-monotonicity of $S_2$ has recently been reported for sheared frictionless 2D ellipses in M. Trulsson, “Rheology and shear jamming of frictional ellipses," J. Fluid Mech. [**849**]{}, 718 (2018). D. B. Nagy, P. Claudin, T. B[ö]{}rzs[ö]{}nyi, and E. Somfai, “Rheology of dense granular flows for elongated particles," Phys. Rev. E [**96**]{}, 062903 (2017).

D. J. Evans and G. P. Morriss, [*Statistical Mechanics of Nonequilibrium Liquids*]{} (ANU Press, 2007).

See Supplemental Material Sec. V for a discussion of detecting overlaps between 2D spherocylinders and 3D ellipsoids, based on the methods of Refs. [@Pournin; @perram1985statistical; @donev2007underconstrained; @donev].

L. Pournin, M. Weber, M. Tsukahara, J.-A. Ferrez, M. Ramaioli, and T. M. Liebling, “Three-dimensional distinct element simulation of spherocylinder crystallization," Granul. Matter [**7**]{}, 119 (2005). J. W. Perram and M. Wertheim, “Statistical mechanics of hard ellipsoids. I. Overlap algorithm and the contact function," Journal of Computational Physics [**58**]{}, 409 (1985). A. Donev, R. Connelly, F. H. Stillinger, and S. Torquato, “Underconstrained jammed packings of nonspherical hard particles: Ellipses and ellipsoids," Phys. Rev. E [**75**]{}, 051304 (2007). A. Donev, “Jammed Packings of Hard Particles," Ph.D. thesis, Princeton University (2006), available from: <http://cims.nyu.edu/~donev/Thesis.pdf>

D. J. Durian, “Foam mechanics at the bubble scale," Phys. Rev. Lett. [**75**]{}, 4780 (1995) and “Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches," Phys. Rev. E [**55**]{}, 1739 (1997). T. A. Marschall, S. V. Franklin, and S. Teitel, “Compression- and shear-driven jamming of U-shaped particles in two dimensions," Granular Matter [**17**]{}, 121 (2015). D. V[å]{}gberg, P. Olsson, and S. Teitel, “Shear banding, discontinuous shear thickening, and rheological phase transitions in athermally sheared frictionless disks," Phys. Rev. E [**95**]{}, 052903 (2017). P. Olsson and S. Teitel, “Critical scaling of shearing rheology at the jamming transition of soft-core frictionless disks," Phys. Rev. E [**83**]{}, 030302(R) (2011). D. V[å]{}gberg, D. Valdez-Balderas, M. A. Moore, P. Olsson, and S. Teitel, “Finite-size scaling at the jamming transition: Corrections to scaling and the correlation-length critical exponent," Phys. Rev. E [****]{} 83, 030303(R) (2011). P. Olsson and S. Teitel, “Herschel-Bulkley shearing rheology near the athermal jamming transition," Phys. Rev. Lett. [**109**]{}, 108001 (2012). P. Olsson, “Dimensionality and viscosity exponent in shear-driven jamming," Phys. Rev. Lett. [**122**]{}, 108003 (2019). See Supplemental Material Sec. I.D for an analysis of the instantaneous $\mathbf{S}_2(\gamma)$ for 2D spherocylinders.

See Supplemental Material Sec. I for tests that our data for 2D spherocylinders with $\alpha=0.001$ are accurate.

For spherocylinders it is ambiguous how to define the location of a contact when two flat sides are touching. Here we take such a contact to have unit strength but distribute its location evenly over the segments of the flat surfaces that are in contact. For small $\alpha$ our results are insensitive to this as compared to other choices.

For 2D spherocylinders we find that near $\phi_J$ the integrated probability for a contact to lie anywhere along the flat sides stays constant as $\alpha\to 0$, $\dot\gamma\to 0$.

See Supplemental Material Sec. III for plots of $\mathcal{P}(\pi/2)$ vs $\phi$ at different $\dot\gamma$.

K. VanderWerf, W. Jin, M. D. Shattuck, and C. S. O’Hern, “Hypostatic jammed packings of frictionless nonspherical particles," Phys. Rev. E [**97**]{}, 012909 (2018). See Supplemental Material Sec. IV for a discussion of the response to a pure shear.

Preliminary results for 2D spherocylinders suggest that our conclusions remain unchanged even with polydispersity in $\alpha$.

**[Orientational Ordering in Athermally Sheared, Aspherical, Frictionless Particles]{}**

**[Supplemental Material]{}**

In this Supplemental Material we provide further details and tests to demonstrate the correctness of our simulations. In Sec. \[valid\] we demonstrate the validity of our results for the nematic order parameter $S_2$ at our smallest asphericity, $\alpha=0.001$, which is key to our conclusion that $S_{2\,\mathrm{max}}$ stays finite as $\alpha\to 0$. In Sec. \[flow\] we discuss translational correlations in our system and demonstrate that there is no smectic ordering into well defined flowing layers. In Sec. \[contact\] we show that the onset for the effect that particle contacts to prefer to lie on the narrowest width of the particles, takes place as the packing $\phi$ increases through the jamming transition. In Sec. \[pure\] we consider the effect of a pure shear deformation on 2D spherocylinders, and contrast with our main results for simple shear. In Sec. \[deltaij\] we provide details on how we determine when two particles are in contact, and compute the corresponding overlap parameters.

Validity of results at $\alpha=0.001$ {#valid}
=====================================

Our argument in the main text, that $\lim_{\alpha\to 0}[S_{2\,\mathrm{max}}]=S_0$ is finite, relied on the assertion that our data at the smallest $\alpha=0.001$ are reliable. In order to argue conversely, i.e., that $S_{2\,\mathrm{max}}$ vanishes as a power law as $\alpha\to 0$, we would have to believe that the value of $S_{2\,\mathrm{max}}$ at $\alpha=0.001$ that is reported in Fig. 4a of the main text is, by some artifact of our simulations, larger than it should be.

Here we provide several tests to support our claim that our data point at $\alpha=0.001$ is indeed correct. Since our simulations for 2D spherocylinders are considerably less time consuming than for 3D ellipsoids, we can make more exacting tests for that case. Hence, here we restrict ourselves to 2D spherocylinders.

Dependence on Shear Strain Rate
-------------------------------

As shown in Fig. 1 of the main text, the nematic order parameter $S_2$ depends on both packing fraction $\phi$ and shear strain rate $\dot\gamma$. However at each $\phi$, $S_2$ approaches a limiting value as $\dot\gamma$ decreases. Here we wish to confirm that we have simulated at small enough $\dot\gamma$ so that the peak value $S_{2\,\mathrm{max}}$ which we find for $\alpha=0.001$ has reached this $\dot\gamma\to 0$ limit. In Fig. \[fSM1\]a we plot $S_2$ vs $\phi$ for our three smallest strain rates $\dot\gamma$, using a system with $N=1024$ particles. Just as was found in Fig. 1 of the main text for a larger $\alpha$, here we see $\dot\gamma$ dependence remaining on the large $\phi$ side of the peak in $S_2$, however there is no $\dot\gamma$ dependence on the low $\phi$ side up to, and including, the peak value. Thus our results of Fig. \[fSM1\]a clearly argue that the value of $S_{2\,\mathrm{max}}$ will not decrease if $\dot\gamma$ were made even smaller.

![Nematic order parameter $S_2$ for $\alpha=0.001$ vs packing fraction $\phi$ for (a) three different shear strain rates $\dot\gamma$ with $N=1024$ particles, and (b) systems with different numbers of particles $N$ for $\dot\gamma=4\times 10^{-7}$. Vertical dashed lines locate the jamming transition of $\alpha=0$ spherical particles, $\phi_J^{(0)}=0.8433$. []{data-label="fSM1"}](fSM1){width="3.5in"}

Dependence on System Size
-------------------------

As one approaches the jamming transition, a correlation length diverges. If one is too close to the jamming transition, finite size effects may become important when the correlation length becomes larger than the length of the system. We thus wish to check that our value of $S_{2\,\mathrm{max}}$ for $\alpha=0.001$ is not affected by such possible finite system size effects. Our results in the main text are for systems with $N=1024$ particles. In Fig. \[fSM1\]b we plot $S_2$ vs $\phi$ at the small strain rate $\dot\gamma=4\times 10^{-7}$, using three different systems sizes with numbers of particles $N=512, 1024$ and 2048. While there is a small dependence on $N$ seen at the largest $\phi$, there is no dependence on $N$ at lower $\phi$ up to and including the peak value. Our value of $S_{2\,\mathrm{max}}$ for $\alpha=0.001$ thus does not suffer from finite size effects.

Dependence on Integration Time Step
-----------------------------------

We should also check if there is any dependence of our results on the size of the finite numerical integration step $\Delta t$. Our results in the main text used a value $\Delta t=0.02 t_0$ with $t_0=D_s^2k_d/k_e$ the unit of time. In Fig. \[fSM2\] we plot $S_2$ vs $\phi$ at the small strain rate $\dot\gamma=4 \times 10^{-7}$, for $\alpha=0.001$, using three different values of the time step $\Delta t/t_0=0.01, 0.02$ and $0.04$. We see that any differences between the data from the three different $\Delta t$ are within the estimated statistical error. We conclude that our integration step of $\Delta t/t_0=0.02$ is small enough to accurately determine $S_{2\,\mathrm{max}}$ for $\alpha=0.001$.

![Nematic order parameter $S_2$ for $\alpha=0.001$ vs packing fraction $\phi$ using different integration time steps $\Delta t/t_0=0.01, 0.02$ and 0.04. The system is sheared at a strain rate $\dot\gamma=4\times 10^{-7}$ and has $N=1024$ particles. Vertical dashed line locates the jamming transition of $\alpha=0$ spherical particles, $\phi_J^{(0)}=0.8433$. []{data-label="fSM2"}](fSM2){width="2.0in"}

Equilibration
-------------

Finally we demonstrate that the rotational degrees of freedom in our system are well equilibrated for our simulations at $\alpha=0.001$. When $\alpha$ is small, the small moment arms of the collisional forces result in small torques, and, depending on the particle density, it can require long shear strains for the rotational degrees of freedom of the system to equilibrate to the proper steady state.

Let us define $S_2(\gamma)$ as the magnitude of the instantaneous nematic order parameter of the individual configuration of the system after it has sheared a strain $\gamma=\dot\gamma t$. For an initial configuration of randomly oriented particles, a system with a finite number of particles $N$ will have some initial value $S_2(0)$. At low densities, where torque inducing collisions are rare, particles will rotate primarily under the influence of the dissipative torque. In this case, since particles with finite $\alpha$ have a non uniform angular velocity that depends on their orientation $\theta_i$, the relative orientations of the particle spines $\boldsymbol{\hat\ell}_i$ will vary with $\gamma$ and hence so will $S_2$. But once the system has strained so that $\gamma=\dot\gamma T$, where $T$ is the period of rotation of an isolated particle, the particles will have returned to their initial orientations and $S_2(\gamma)$ will have returned to its initial value $S_2(0)$. We thus expect to see an oscillating $S_2(\gamma)$ with period $\dot\gamma T$. We have observed such behavior for small $\alpha$ at low densities. However, as the density increases the rate of collisions increases. These collisions will perturb this oscillatory behavior until, after a sufficiently large strain has been applied, the particle orientations have lost memory of their initial configuration. The particle orientations will then sample a stationary steady state distribution. $S_2(\gamma)$ will then be constant, aside from random fluctuations that we might expect should decrease as $1/\sqrt{N}$ as the number of particles $N$ increases.

In Fig. \[fSM3\]a we plot $S_2(\gamma)$ vs $\gamma$ for a shear strain rate $\dot\gamma=4\times 10^{-7}$ at a packing $\phi=0.838$ near the peak in $S_2$, for a system with $N=1024$ particles with $\alpha=0.001$. We see that $S_2(\gamma)$ appears, as desired, to consist of random fluctuations about a constant average. The dashed horizontal line in Fig. \[fSM3\]a is the average $\langle S_2(\gamma)\rangle=(1/\Delta\gamma)\int_{\gamma_i}^{\gamma_f}d\gamma\,S_2(\gamma)$, where $\Delta\gamma=\gamma_f-\gamma_i$; we start averaging only after an initial shear strain of $\gamma_i=10$ so as to avoid any initial transients, and average up to a final $\gamma_f=150$. The solid horizontal line represents the ensemble average $S_2$, as considered elsewhere in this work. To be clear, $S_2(\gamma)$ is the largest eigenvalue of the orientational ordering tensor $T_{\mu\nu}(\gamma)$ as computed for the individual configuration at strain $\gamma$, while $S_2$ is the largest eigenvalue of the orientational ordering tensor $\langle T_{\mu\nu}\rangle$ averaged over the length of the shearing run from $\gamma_1$ to $\gamma_2$. Since the relation between eigenvalue and tensor is not linear, these two averages of $S_2$ need not be equal, and in Fig. \[fSM3\]a we see that there is indeed a small difference. Since the direction of the nematic director is optimized to give the largest possible $S_2$, and since the direction of the nematic director obtained from $T_{\mu\nu}(\gamma)$ fluctuates as $\gamma$ varies from configuration to configuration (as opposed to the director obtained from $\langle T_{\mu\nu}\rangle$ which is fixed), we expect that $\langle S_2(\gamma)\rangle$ will be somewhat larger than $S_2$, and this is indeed what is observed in Fig. \[fSM3\]a. In Fig. \[fSM3\]b we plot $\langle S_2(\gamma)\rangle - S_2$ vs $N$ and see that this difference is going to zero as $N$ increases. In the same figure we also plot the standard deviation $\sigma_{S_2(\gamma)}= \sqrt{\langle S_2^2(\gamma)\rangle -\langle S_2(\gamma)\rangle^2}$ vs $N$ and see that it also vanishes as $N$ increases.

![(a) Instantaneous nematic order parameter $S_2(\gamma)$ vs shear strain $\gamma$, for $\alpha=0.001$ and shear strain rate $\dot\gamma=4\times 10^{-7}$ at packing fraction $\phi=0.838$ near the peak value $S_{2\,\mathrm{max}}$. The horizontal dashed line is the average over these instantaneous values $\langle S_2(\gamma)\rangle$, while the horizontal solid line is $S_2$ as obtained from averaging the orientational ordering tensor over the entire shearing run. The system has $N=1024$ particles. (b) Difference $\langle S_2(\gamma)\rangle - S_2$ vs number of particles $N$, and standard deviation $\sigma_{S_2(\gamma)}$ vs $N$; the dashed line is $\sim 1/\sqrt{N}$ for comparison. []{data-label="fSM3"}](fSM3){width="3.5in"}

Next, we consider the Fourier transform of $S_2(\gamma)$ in order to check that the frequency spectrum of the fluctuating noise seen in Fig. \[fSM3\]a is broad without any peaks that could indicate vestigial oscillations due to poor equilibration. Since $S_2(\gamma)$ is plotted in terms of the dimensionless time $\gamma=\dot\gamma t$, in Fig. \[fSM4\] we plot the Fourier transform $\mathcal{F}[S_2]$ as a function of the dimensionless frequency $\omega/\dot\gamma$. We see that the spectrum is indeed broad with no peaks. The high frequency tail is roughly power law with an exponent $~1.3$, however that exponent changes a bit depending on the range of $\gamma$ that is used in the fit.

![Fourier transform of $S_2(\gamma)$, $\mathcal{F}[S_2]$, vs dimensionless frequency $\omega/\dot\gamma$. The high frequency tail is fit to an inverse power law (dashed line) and gives an exponent $\sim 1.3$. []{data-label="fSM4"}](fSM4){height="1.75in"}

![(a) Instantaneous nematic director angle $\theta_2(\gamma)$ vs shear strain $\gamma$, for $\alpha=0.001$ and shear strain rate $\dot\gamma=4\times 10^{-7}$ at packing fraction $\phi=0.838$ near the peak value $S_{2\,\mathrm{max}}$. The horizontal dashed line is the average over these instantaneous values $\langle \theta_2(\gamma)\rangle$, while the horizontal solid line is $\theta_2$ as obtained from averaging the orientational ordering tensor over the entire shearing run. The system has $N=1024$ particles. (b) Fourier transform of $\theta_2(\gamma)$, $\mathcal{F}[\theta_2]$, vs dimensionless frequency $\omega/\dot\gamma$. The high frequency tail is fit to an inverse power law (dashed line) and gives an exponent $\sim 1.5$. []{data-label="fSM5"}](fSM5){width="3.5in"}

Lastly we consider a similar analysis of the orientation angle $\theta_2$ of the nematic director. In Fig. \[fSM5\]a we plot the instantaneous $\theta_2(\gamma)$ vs $\gamma$ for the same parameters as in Fig. \[fSM3\]a, $\dot\gamma=4\times 10^{-7}$, $\phi=0.838$, $\alpha=0.001$. We see what appear to be random fluctuations about a constant average value. The dashed horizontal line is the average $\langle \theta_2(\gamma)\rangle = (1/\Delta\gamma)\int_{\gamma_i}^{\gamma_f}d\gamma\,\theta_2(\gamma)$, while the solid horizontal line is $\theta_2$ obtained from the ensemble averaged orientation tensor $\langle T_{\mu\nu}\rangle$. In Fig. \[fSM5\]b we plot the Fourier transform $\mathcal{F}[\theta_2]$ vs the dimensionless frequency $\omega/\dot\gamma$. We see a broad spectrum with a power law tail decreasing with an exponent $\sim 1.5$ (the exact value of this exponent is sensitive to the range of data used in the fit). There are no peaks in $\mathcal{F}[\theta_2]$ to indicate any oscillatory motion, thus giving support to the assertion in the main text that, while individual particles tumble with an average angular velocity $\langle \omega_i\rangle$, there is no coherent tumbling of the nematic order parameter $\mathbf{S}_2$. Our results in this section thus confirm that our spherocylinder simulations at $\alpha=0.001$ are indeed well equilibrated.

Spatial Correlations {#flow}
====================

Since our system has finite nematic orientational order, we wish to check whether there might also be smectic translational order, with particles flowing in well defined layers oriented in the direction of the flow. To test for this, we measure the following transverse correlation function of the particle center of mass density $n(\mathbf{r})$. We first define, $$n(y)=\dfrac{1}{\Delta y \mathcal{L}_\perp^{d-1}}\int_{y-\Delta y/2}^{y+\Delta y/2}\!\!\! dy^\prime\int_0^{\mathcal{L}_\perp}\!\!\!d\mathbf{r}_\perp \, n(y^\prime, \mathbf{r}_\perp).$$ $n(y)$ is just he number of particles per unit volume whose center of mass lies in a layer of small width $\Delta y$ that spans the system in the orthogonal directions. For $d=2$ dimensions, $\mathbf{r}_\perp = x\mathbf{\hat x}$ and $\mathcal{L}_\perp=\mathcal{L}_x$, the length of the system in the $\mathbf{\hat x}$ direction; for $d=3$, $\mathbf{r}_\perp=x\mathbf{\hat x}+z\mathbf{\hat z}$ and $\mathcal{L}_\perp=\mathcal{L}_x=\mathcal{L}_z$. We then define the correlation $$C(y) = \dfrac{\mathcal{L}_\perp^{d-1}}{n\mathcal{L}_y}\int_0^{\mathcal{L}_y}\! \!dy^\prime \left[\langle n(y+y^\prime)n(y^\prime)\rangle - \langle n\rangle^2\right],$$ where the prefactor is chosen so that $C(y)$ is independent of the system size.

![Spatial correlations in the direction of the flow velocity gradient, $C(y)$ vs $y$, for 2D spherocylinders. (a) Results for small $\alpha$, including $\alpha=0$, at the common value of $\phi=0.845$, just above the jamming fraction for circular disks, $\phi_J^{(0)}=0.8433$; strain rates $\dot\gamma$ are as indicated in the figure. (b) Results for larger values of $\alpha$, at low strain rates $\dot\gamma$, just above their respective jamming fractions $\phi_J(\alpha)$; values of $\dot\gamma$ and $\phi$ are indicated in the figure. []{data-label="C-vs-y-2D"}](C-vs-y-2D){width="3.5in"}

We consider first the case of spherocylinders in 2D, where we average over large total strains $\gamma\approx 130$, thus allowing for accurate measurements of $C(y)$. For our calculations we use a layer width $\Delta y = 0.01$ for $\alpha\le 0.01$, and $\Delta y=0.025$ for larger $\alpha$. In Fig. \[C-vs-y-2D\]a we plot $C(y)$ vs $y$ for $\alpha=0$, $0.001$, and $0.01$ at $\phi=0.845$, which is just slightly above the jamming $\phi_J^{(0)}=0.8433$ for circular disks. In each case we use the smallest $\dot\gamma$ we have simulated at each $\alpha$, i.e. $10^{-6}$, $4\times10^{-7}$ and $4\times 10^{-7}$ respectively. We see that the $C(y)$ for these three cases are almost indistinguishable; there is nothing that signals a singular behavior as $\alpha\to 0$. We see sharp peaks at $y=1$, 1.2, and 1.4, which are the nearest neighbor separations for just contacting small-small, small-big, and big-big pairs. At larger $y$ we see oscillations with a period of 1.2, the average spacing between contacting particles. However these oscillations clearly decay to zero as $y$ increases, thus demonstrating that there is only short ranged order in the direction of the flow velocity gradient. Fitting the heights of the larger $y$ peaks to an exponential, we find a decay length between $1$ and $2$.

In Fig. \[C-vs-y-2D\]b we plot $C(y)$ vs $y$ for larger values of $\alpha$, at our lowest strain rate for each case, and at a packing fraction $\phi$ that is slightly above the respective jamming fraction $\phi_J$ for each $\alpha$. We again see similar behavior: oscillations that decay to zero as $y$ increases. As $\alpha$ increases, and the particles become increasingly non-spherical, the sharp peaks near $y=1$, 1.2 and 1.4 broaden and the peaks at $y>2$ shift to slightly larger values of $y$; results for $\alpha=1$ and $\alpha=4$ are nearly indistinguishable for $y>2$. However the average spacing between peaks remains $\sim 1.2$ and the decay length remains in the range 1 to 2. We have verified that similar behavior occurs as either $\phi$ or $\dot\gamma$ is varied. We thus conclude that particles do not flow in well defined, spatially ordered, layers and so there is no smectic ordering.

Our 3D simulations are much more time consuming and we only shear to total strains $\gamma\approx 1.4$, thus greatly reducing the number of independent samples we have to average over when computing $C(y)$. To keep statistical accuracy reasonable, we therefore average over a thicker (as compared to 2D) layer of width $\Delta y = 0.18$ to define $n(y)$, so as to have more particles in the layer and so smaller fluctuations. Our results for the correlations of 3D ellipsoids are shown in Fig. \[C-vs-y-3D\]. While the larger $\Delta y$ means we lack the finer scale features seen in Fig. \[C-vs-y-2D\] for 2D, we continue to see similar decaying oscillations, characteristic of the absence of any long range translational ordering.

![Spatial correlations in the direction of the flow velocity gradient, $C(y)$ vs $y$, for 3D ellipsoids. (a) Results for small $\alpha$, including $\alpha=0$, at strain rate $\dot\gamma=10^{-7}$ and the common value of $\phi=0.654$, just above the jamming fraction for spheres, $\phi_J^{(0)}=0.649$. (b) Results for larger values of $\alpha$, at strain rates $\dot\gamma=5\times 10^{-7}$, just above their respective jamming fractions $\phi_J(\alpha)$; values of $\phi$ are indicated in the figure. []{data-label="C-vs-y-3D"}](C-vs-y-3D){width="3.5in"}

Note, since we use $N=1024$ particles in both 2D and 3D, the system length for our 2D systems is $\mathcal{L}\sim 37$, while for 3D it is $\mathcal{L}\sim 11$. Thus in 3D the oscillations in $C(y)$ have not quite decayed to zero before one reaches $y=\mathcal{L}/2$, where the periodic boundaries influence our results and give a larger $C(y)$ than would be found in a larger system. Nevertheless our results in 3D are consistent with decaying correlations, and so the absence of any smectic ordering.

Contact location distribution {#contact}
=============================

In the main text we showed in Fig. 5 that the probability density per unit surface area $\mathcal{P}(\vartheta)$, for a particle to have a contact at polar angle $\vartheta$ on its surface, had a sharp peak at $\vartheta=\pi/2$, where the particle width is narrowest. The height of this peak increases as the asphericity $\alpha$ decreases. The results for $\mathcal{P}(\vartheta)$ vs $\vartheta$ shown in Fig. 5 were for a small strain rate $\dot\gamma$ at a fixed packing fraction near the jamming transition for spherical particles, $\phi\approx \phi_J^{(0)}$.

In Fig. \[fSM6\] we plot the peak height $\mathcal{P}(\pi/2)$ vs packing $\phi$ at fixed small $\alpha$, for different values of $\dot\gamma$. In (a) we show 2D spherocylinders at $\alpha=0.03$ and in (b) 3D ellipsoids at $\alpha=0.05$. We see that as $\dot\gamma$ decreases, $\mathcal{P}(\pi/2)$ increases to a limiting curve, which rises rapidly as $\phi$ approaches $\phi_J^{(0)}$, and then stays above the spherical particle value of unity as $\phi$ increases above the jamming transition. Thus the onset for the contacts to preferentially lie along the narrowest width of the particle takes place as $\phi$ passes through the jamming transition.

![Peak probability $\mathcal{P}(\pi/2)$ vs packing $\phi$ for different strain rates $\dot\gamma$ for (a) 2D spherocylinders at $\alpha=0.03$ and (b) 3D ellipsoids at $\alpha=0.05$. As $\dot\gamma$ decreases, the peak value $\mathcal{P}(\pi/2)$ increases until it saturates. Vertical dashed lines denote the jamming point of spherical particles $\phi_J^{(0)}$, while horizontal dashed lines indicate the value of unity expected for a spherical particle. []{data-label="fSM6"}](fSM6){width="3.5in"}

Response to a pure shear deformation {#pure}
====================================

It is interesting to compare the response of our system to a pure, rather than a simple, shear. For simplicity we consider this for the case of our 2D spherocylinders. In this model the dissipative force is determined by the relative velocity of the particle with respect to an affinely deformed background host medium. We define the local velocity $\mathbf{v}_\mathrm{host}(\mathbf{r})$ of this background host in terms of a strain rate tensor $\dot{\boldsymbol{\Gamma}}$, i.e., $\mathbf{v}_\mathrm{host}=\dot{\boldsymbol{\Gamma}}\cdot\mathbf{r}$. A simple shear deformation can be decomposed into the sum of a pure shear and a uniform rotation, $$\dot{\boldsymbol{\Gamma}}_\mathrm{ss}=\dot{\boldsymbol{\Gamma}}_\mathrm{ps}+\dot{\boldsymbol{\Gamma}}_\mathrm{rot}.$$ For our coordinate system with simple shear flow in the $\mathbf{\hat x}$ direction, this becomes, $$\left[
\begin{array}{cc}
0 & \dot\gamma\\
0 & 0
\end{array}
\right]
=
\left[
\begin{array}{cc}
0 & \dot\gamma/2\\
\dot\gamma/2 & 0
\end{array}
\right]
+\left[
\begin{array}{cc}
0 & \dot\gamma/2\\
-\dot\gamma/2 & 0
\end{array}
\right].$$ The first term on the right hand side is a pure shear, with expansion along the $(1,1)$ diagonal and compression along the $(1,-1)$ diagonal, both at the rate $\dot\gamma/2$ so as to keep the area constant. The second term is a clockwise rotation $(-\dot\gamma/2)\mathbf{\hat z}\times \mathbf{r}$, with angular velocity $-\dot\gamma/2$. It is this second term which drives the continuous rotation of particles under simple shear, resulting in the finite $-\langle\omega_{zi}\rangle/\dot\gamma >0$ seen in Fig. 3 of the main text.

Under pure shear there is no such rotational drive, and particles try to relax from their initial orientation to one aligned with the expansive direction of the pure shear. Rotating coordinates so that the expansive direction is $\mathbf{\hat x}$ and the compressive direction is $\mathbf{\hat y}$, the rotational equation of motion for pure shear becomes, $$\dot\theta_i=-(\dot\gamma/2)[\Delta I_i/I_i]\sin 2\theta_i + \tau_i^\mathrm{el}/(k_dv_iI_i).$$ For an isolated particle where $\tau_i^\mathrm{el}=0$, particles will exponentially relax to $\theta_i=0$ or $\pi$ with a relaxation time $t_0$ set by the total strain $\gamma_0=\dot\gamma t_0 = I_i/\Delta I_i$. Thus, at low $\phi$ near this isolated particle limit, we expect to find near perfect nematic ordering with $S_2\approx 1$ and $\theta_2=0$. However, as the asphericity $\alpha$ of the particles vanishes, the relaxation time needed to achieve this highly ordered state diverges as $\dot\gamma t_0=(I_i/\Delta I_i)\sim 1/\alpha$.

To investigate the response to pure shear at dense $\phi$, we have carried out numerical simulations. A practical limitation of pure shear simulations is that, unlike for simple shear, there is a limit to the total strain $\gamma$ that can be applied to a finite numerical system before the system collapses to a narrow height of order one particle width, $L_y(\gamma)=L_y(0)\mathrm{e}^{-\gamma/2}\sim O(1)$. To increase the total possible strain $\gamma$, we use systems of $N=1024$ particles with an initial system aspect ratio of $L_y(0)/L_x(0)=8$, and shear to a strain $\gamma$ such that $L_y(\gamma)/L_x(\gamma)=1/8$, thus allowing a maximum strain of $\gamma=\ln 64\approx 4.2$. We use a strain rate $\dot\gamma=10^{-6}$, and average over four independently generated samples.

![(a) Magnitude $S_2$ and (b) direction $\theta_2$ of the nematic order parameter vs pure strain $\gamma=\dot\gamma t$, at different packing fractions $\phi$ for nearly circular particles with $\alpha=0.001$ at strain rate $\dot\gamma=10^{-6}$. A sparse set of symbols is used to help differentiate curves of different $\phi$, with many data points existing between adjacent symbols on any curve. Representative error bars are shown at integer values of $\gamma$. []{data-label="S2-vs-g"}](S2-vs-g){width="3.5in"}

Here we present results for nearly circular particles at our smallest $\alpha=0.001$. In Fig. \[S2-vs-g\] we plot the nematic order parameter magnitude $S_2$ and orientation $\theta_2$ vs pure shear strain $\gamma$, for several different packing fractions $\phi$. As $\gamma$ increases, $S_2$ increases and $\theta_2$ decays from its initial random value to zero, in agreement with the expectation that particles try to relax to their preferred orientation aligned with the expansive direction $\mathbf{\hat x}$. However we see that we are only able to reach the desired steady state, where $S_2$ plateaus to a constant value as $\gamma$ increases, for relatively dense systems close to and above jamming, $\phi \ge 0.84$.

![Magnitude of the nematic order parameter $S_2$ vs packing $\phi$, comparing pure shear with simple shear, for nearly circular particles with $\alpha=0.001$. For pure shear the strain rate is $\dot\gamma=10^{-6}$, while for simple shear $\dot\gamma=10^{-7}$. Results represent steady state values, except for the pure shear case at $\phi=0.835$ where steady-state has not quite been reached; the value shown at this $\phi$ is therefore a lower bound on the steady state limit. []{data-label="S2-vs-phi"}](S2-vs-phi){width="3.2in"}

In Fig. \[S2-vs-phi\] we show the resulting steady state values of $S_2$ vs packing $\phi$, comparing results from simple shear with those from pure shear. We see dramatically different behavior at low $\phi$. For simple shear the nearly uniform rotation of the $\alpha=0.001$ particles results in a small $S_2$, while for pure shear the relaxation to the expansive direction gives a large $S_2$. As $\phi$ increases, so does the rate of particle collisions. For pure shear the collisions and resulting excluded volume inhibit perfect alignment of particles and $S_2$ decreases. For simple shear the increasing collisions initially cause the rotation to slow (see Fig. 3a of the main text) and consequently $S_2$ to increase, but upon further increasing $\phi$ towards $\phi_J$ and going above, excluded volume effects similar to that in pure shear presumably inhibit alignment and cause $S_2$ to decrease, and we find that $S_2$ for both pure and simple shear become comparable and behave similarly. The non-monotonic behavior of $S_2$ in simple shear is thus a consequence of the rotational drive, present in simple shear but absent in pure shear. However in both simple and pure shear, we find that $S_2$ at jamming remains surprisingly large, even though the particles are extremely close to circular, with the flat sides of the spherocylinders comprising only a fraction $\alpha/(\alpha+\pi/2)=6.4\times 10^{-4}$ of the particle perimeter.

Determination of contacts and overlaps {#deltaij}
======================================

In this section we summarize how we determine if two particles are overlapping, and if so, what is the point of contact between them. For our 2D spherocylinders, we use the efficient algorithm of Pournin et al. [@Pournin] to compute the shortest distance $r_{ij}$ between the spines of two spherocylinders $i$ and $j$. The line of length $r_{ij}$ that connects the two spines we will call the line IJ. Whenever $r_{ij}<d_{ij}=(D_i+D_j)/2$, with $D_i$ the diameter of the endcap of spherocylinder $i$, the two spherocylinders are overlapping. We then define the point of contact ${\mathbf{r}}_C$, at which the elastic force acts, as the distance $[D_i/(D_i+D_j)]r_{ij}$ from the spine of spherocylinder $i$, along the line IJ.

For our 3D ellipsoids, the procedure is more complicated. As illustrated in Fig. \[overlap\], for two overlapping ellipsoids $i$ and $j$ one can define a scale factor $\delta_{ij} < 1$ such that there exists a unique point of contact ${\mathbf{r}}_C$ between these ellipsoids when their axes are rescaled by the common factor $\delta_{ij}$, keeping their center of mass positions fixed. This scale factor $\delta_{ij}$ can be computed using a method introduced by Perram and Wertheim (PW) [@perram1985statistical] which has been applied to the study of jammed packings of ellipsoidal particles [@donev2007underconstrained; @donev]. Here we briefly summarize this method.

For any position ${\mathbf{r}}$, we define the scale function $\delta_i({\mathbf{r}})$ such that ${\mathbf{r}}$ will lie on the surface of ellipsoid $i$ if its axes are rescaled by $\delta_i({\mathbf{r}})$. We then introduce the contact function $F({\mathbf{r}},\lambda)$ defined for two ellipsoids $i$ and $j$, $$F({\mathbf{r}}, \lambda) = \lambda \delta_i^2({\mathbf{r}}) + (1 - \lambda) \delta_j^2({\mathbf{r}}),
  \label{eq:contact_function}$$ where $\lambda\in [0,1]$. It has then been demonstrated [@perram1985statistical] that there exists an ${\mathbf{r}}(\lambda)$ such that $$\label{eq:nablaF}
  {\boldsymbol{\nabla}}F({\mathbf{r}}(\lambda), \lambda) = 0,$$ where ${\boldsymbol{\nabla}}\equiv \partial/\partial\boldsymbol{r}$. This implies that $$\lambda {\boldsymbol{\nabla}}\delta_i^2({\mathbf{r}}(\lambda)) = -(1 - \lambda) {\boldsymbol{\nabla}}\delta_j^2({\mathbf{r}}(\lambda)),
  \label{eq:nabladeltai}$$ which shows that when ellipsoids $i$ and $j$ are rescaled by factors $\delta_i({\mathbf{r}}(\lambda))$ and $\delta_j({\mathbf{r}}(\lambda))$ respectively, the point ${\mathbf{r}}(\lambda)$ lies on the surfaces of both ellipsoids, and the normal vectors to the surfaces at this point are parallel but pointing in opposite directions, so that the two ellipsoids are tangent at ${\mathbf{r}}(\lambda)$.

PW further showed [@perram1985statistical] that $F({\mathbf{r}}(\lambda), \lambda)$, as a function of $\lambda\in[0,1]$, has a unique maximum at $\lambda^*$, such that $$\begin{aligned}
       0 & =\left. \frac{d F({\mathbf{r}}(\lambda), \lambda)}{d\lambda}  \right|_{\lambda=\lambda^*}\\[10pt]
    & =\left. \frac{\partial F({\mathbf{r}}(\lambda), \lambda)}{\partial \lambda}\right|_{\lambda=\lambda^*}  \!\!+ {\mathbf{r}}(\lambda^*) \cdot {\boldsymbol{\nabla}}F({\mathbf{r}}(\lambda^*), \lambda^*),
  \end{aligned}$$ where the second term vanishes due to [Eq. (\[eq:nablaF\])]{}. From [Eq. (\[eq:contact\_function\])]{} we then find $\delta_i^2({\mathbf{r}}(\lambda^*)) = \delta_j^2({\mathbf{r}}(\lambda^*))$, which means that the scale factor is the same for both ellipsoids, and $$\label{eq:hmm}
  \delta_i^2({\mathbf{r}}(\lambda^*)) = \delta_j^2({\mathbf{r}}(\lambda^*)) = F({\mathbf{r}}(\lambda^*), \lambda^*).$$ The scale factor $\delta_{ij}$ that we are seeking is thus defined as $$\delta_{ij}^2 = \max_{ \lambda \in [0,1]} [ F({\mathbf{r}}(\lambda), \lambda)].
  \label{eq:scale_factor_max}$$

With this notation, we define the point of contact between ellipsoids $i$ and $j$ as ${\mathbf{r}}_C = {\mathbf{r}}(\lambda^*)$. It is thus the unique point common to ellipsoids $i$ and $j$ when both are rescaled with a common factor $\delta_{ij}$.

(0,0) circle (0.5pt); (0,-0.03) node\[scale=1.5, below\][$\mathbf{r}_i$]{}; (0,0) ellipse (1 and 0.5); (0.8,0.4) circle (0.5pt); (0.8,0.4+0.02) node\[scale=1.5, above\][$\mathbf{r}_j$]{}; (0.8,0.4) ellipse (1 and 0.5); (0,0) ellipse (0.655762\*1 and 0.655762\*0.5); (0.8,0.4) ellipse (0.655762\*1 and 0.655762\*0.5); (0.555670,0.174160) circle (0.5pt) ; (0.555670+0.01,0.174160+0.01) node\[scale=1.5, right\][$\mathbf{r}_C$]{};

To compute the scale factor $\delta_{ij}$ defined in [Eq. (\[eq:scale\_factor\_max\])]{}, we use a method derived from Ref. [@donev]. An ellipsoid $i$ is defined by its center of mass position ${\mathbf{r}}_i$, the lengths of its axes $(a_1, a_2, a_3)$, and the rotation matrix $\boldsymbol{\mathcal{Q}}_i$ that rotates the $(x, y, z)$ directions of the lab coordinate frame onto the principal axes of the ellispoid. We then introduce the matrix, $${\mathbf{B}}_i = \boldsymbol{\mathcal{Q}}_i\cdot
  \begin{pmatrix} a^{-2}_1 & 0 & 0 \\ 0 & a^{-2}_2 & 0 \\ 0 & 0 & a^{-2}_3 \end{pmatrix}\cdot
  \boldsymbol{\mathcal{Q}}_i^{-1},$$ which is symmetric due to the orthogonal nature of $\boldsymbol{\mathcal{Q}}_i$, and gives an explicit definition of the scale function $\delta_i({\mathbf{r}})$, $$\delta_i^2({\mathbf{r}}) = ({\mathbf{r}}- {\mathbf{r}}_i)\cdot {\mathbf{B}}_i\cdot ({\mathbf{r}}- {\mathbf{r}}_i).
  \label{eq:B}$$ [Eq. (\[eq:nabladeltai\])]{} then becomes $$\lambda {\mathbf{B}}_i\cdot({\mathbf{r}}(\lambda) - {\mathbf{r}}_i) = - (1 - \lambda) {\mathbf{B}}_j\cdot({\mathbf{r}}(\lambda) - {\mathbf{r}}_j).
  \label{eq:BiBj}$$ After introducing [@perram1985statistical] $${\mathbf{Y}}_{ij}(\lambda) = \lambda {\mathbf{B}}_j^{-1} + (1 - \lambda) {\mathbf{B}}_i^{-1},$$ and defining ${\mathbf{r}}_{ji}={\mathbf{r}}_j-{\mathbf{r}}_i$, [Eq. (\[eq:BiBj\])]{} gives expressions for the distances between the contact point ${\mathbf{r}}(\lambda)$ and the centers of the ellipsoids, $$\begin{aligned}
    {\mathbf{r}}(\lambda) - {\mathbf{r}}_i &= (1 - \lambda) {\mathbf{B}}_i^{-1}\cdot{\mathbf{Y}}^{-1}_{ij}(\lambda)\cdot{{\mathbf{r}}_{ji}},\\
    {\mathbf{r}}(\lambda) - {\mathbf{r}}_j &= -\lambda {\mathbf{B}}_j^{-1}\cdot{\mathbf{Y}}^{-1}_{ij}(\lambda)\cdot{{\mathbf{r}}_{ji}}.
  \end{aligned}
  \label{eq:rlambda-r}$$ As discussed above, the unique contact point ${\mathbf{r}}_C$ for equal scale factors, $\delta_{ij}=\delta_i({\mathbf{r}}_C)=\delta_j({\mathbf{r}}_C)$, is found by maximizing the contact function $F({\mathbf{r}}(\lambda),\lambda)$ with respect to $\lambda$. Using the above results and [Eq. (\[eq:B\])]{} in [Eq. (\[eq:contact\_function\])]{} thus gives, $$\begin{aligned}
\nonumber
  F({\mathbf{r}}(\lambda), \lambda) & = &
  \lambda (1 - \lambda) {{\mathbf{r}}_{ji}}\cdot {\mathbf{Y}}^{-1}_{ij}(\lambda)\cdot {{\mathbf{r}}_{ji}}\qquad\qquad
  \\[10pt]
  & = & \frac{\lambda(1 - \lambda){{\mathbf{r}}_{ji}}\cdot \mathrm{adj}[{\mathbf{Y}}_{ij}(\lambda)]
    \cdot{{\mathbf{r}}_{ji}}}{\mathrm{det}[{\mathbf{Y}}_{ij}(\lambda)]}\\[8pt]\nonumber
  &\equiv& \frac{p_{ij}(\lambda)}{q_{ij}(\lambda)},\end{aligned}$$ where $\mathrm{adj}[\ldots]$ denotes the adjugate matrix (whose element $(\alpha,\beta)$ is equal to the determinant of the $2\times2$ sub-matrix obtained after eliminating row $\beta$ and column $\alpha$ from the original $3\times 3$ matrix), and $\mathrm{det}[\ldots]$ denotes the determinant. The functions $p_{ij}(\lambda)$ and $q_{ij}(\lambda)$ are polynomials in $\lambda$ of degree 4 and 3 respectively.

The desired parameter $\lambda^*$, at which $F({\mathbf{r}}(\lambda), \lambda)$ is maximized, is then the unique root in the interval $[0, 1]$ of the 6th degree polynomial $$h_{ij}(\lambda) = p'_{ij}(\lambda)q_{ij}(\lambda) - p_{ij}(\lambda)q'_{ij}(\lambda),$$ i.e., $h_{ij}(\lambda^*)=0$, where primes above denote derivatives with respect to $\lambda$.

Finally, to determine ellipsoid elastic interactions, we investigate all pairs of ellipsoids whose center of mass separation $|{\mathbf{r}}_i-{\mathbf{r}}_j|$ is small enough that the ellipsoids might be overlapping. We then apply the above procedure to determine $\delta_{ij}$. If the resulting $\delta_{ij}>1$, then the pair of ellipsoids are in fact not overlapping and so have no interaction. If $\delta_{ij}\le 1$, then the ellipsoids overlap and the point of contact is taken as ${\mathbf{r}}_C$.