Two-temperature activity induces liquid-crystal phases inaccessible in equilibrium

The equilibrium liquid-crystalline properties of anisotropic particles are well understood Onsager (1949); De Gennes and Prost (1993); Bolhuis and Frenkel (1997); McGrother et al. (1996); Cuetos et al. (2002); Cuetos and Martínez-Haya (2015); Maiti et al. (2002); Lansac et al. (2003); Thirumalai (1986). Onsager Onsager (1949) first developed an analytical theory to describe the transition from the isotropic (I) phase to the nematic (N) phase, which has uniaxial apolar orientational order, and predicted that a phase with purely nematic order, breaking no other symmetries, cannot arise for hard rods with aspect ratio $L/D<3.70$ Onsager (1949); Bolhuis and Frenkel (1997); De Gennes and Prost (1993). In this Letter, we inquire into the extension of Onsager’s limit to active matter, in the specific context of two-temperature systems.

Active matter is driven locally by a constant supply of free energy to its constituent particles, which dissipate it by performing mechanical work Marchetti et al. (2013); Ramaswamy (2010); Toner et al. (2005); Romanczuk et al. (2012); Prost et al. (2015); Ramaswamy (2019); Gompper et al. (2020); Aditi Simha and Ramaswamy (2002); Bowick et al. (2022); Wensink et al. (2012); Doostmohammadi et al. (2018); Thijssen et al. (2020); Maitra and Voituriez (2020); Nejad and Yeomans (2022); Buttinoni et al. (2013); Bialké et al. (2015); Ivlev et al. (2015); Bechinger et al. (2016); Mandal et al. (2019); Löwen (2020); Speck et al. (2014); McCandlish et al. (2012); Redner et al. (2013). In flocking models, arguably the most familiar examples, activity is linked to a vector order parameter Marchetti et al. (2013); Vicsek and Zafeiris (2012); Ramaswamy (2010); Toner and Tu (1995); Toner et al. (2005). In scalar active matterWittkowski et al. (2014); Saha et al. (2020); Grosberg and Joanny (2015); Weber et al. (2016); Ivlev et al. (2015), activity enters by minimally breaking detailed balance in scalar Halperin-Hohenberg models Hohenberg and Halperin (1977) or, as in our present work, by introducing two (or more) species of particles coupled to thermal baths at distinct temperatures Grosberg and Joanny (2015, 2018); Ilker and Joanny (2020); Weber et al. (2016); Ganai et al. (2014); Smrek and Kremer (2017, 2018); Chari et al. (2019); Chattopadhyay et al. (2021); V et al. (2021). The temperatures in question should be viewed not as thermodynamic but emergent, arising from the effective diffusivities of the multiple motile species. Two-temperature models have accounted for chromatin organization in the cell nucleus Ganai et al. (2014); phase-separation and self-organized structure formation in binary mixtures of Brownian soft disks Weber et al. (2016) and Lennard-Jones (LJ) particles Chari et al. (2019), and in polymer systems Smrek and Kremer (2017, 2018). In a recent work Chattopadhyay et al. (2021), we have implemented this idea in a system of soft repulsive spherocylinders (SRSs) of aspect ratio $L/D=5$ (where $L$ and $D$ are the effective length and diameter defined by the anisotropic repulsive potential Allen et al. (1993); Vega and Lago (1994); Bolhuis and Frenkel (1997); Cuetos et al. (2002); Earl et al. (2001)) and showed that increasing the temperature of the hot particles promotes liquid-crystal ordering in the cold particles, shifting the I-N phase boundary to lower densities than its equilibrium location. Here, we aim to explore ordering transitions of SRS of different $L/D$, in particular, those below Onsager’s limit.

To do so, we have carried out a series of Molecular Dynamics (MD) simulations of a system of SRS with $L/D=5,3$ and $2$. The model and simulation protocol are the same as those in our previous work Chattopadhyay et al. (2021). Here we give brief details for completeness. The SRS interact through the Weeks-Chandler-Andersen potential Weeks et al. (1971):

$$\begin{split}U_{SRS}&=4\varepsilon\left[\left(\dfrac{D}{d_{m}}\right)^{12}-\left(\dfrac{D}{d_{m}}\right)^{6}\right]+\epsilon\qquad\hbox{if}\quad d_{m}<2^{\frac{1}{6}}D\\ &=0\qquad\qquad\qquad\qquad\qquad\qquad\qquad\hbox{if}\quad d_{m}\geq 2^{\frac{1}{6}}D.\end{split}$$

(1)

Here $d_{m}$, the shortest distance between two spherocylinders, implicitly determines their relative orientation and the direction of the interaction force Allen et al. (1993); Vega and Lago (1994); Bolhuis and Frenkel (1997); Cuetos et al. (2002); Earl et al. (2001). We build the initial configuration in a hexagonal close-packed (HCP) crystal structure and perform MD simulations at constant particle number, pressure and temperature (NPT) with periodic boundary conditions in all three directions Rotunno et al. (2004); Maiti et al. (2002); Lansac et al. (2003). For each aspect ratio, we simulate a wide range of pressures spanning the transition from the crystal to the isotropic phase and characterize the phases by calculating the nematic order parameter and appropriate pair correlation functions. For a system of $N$ spherocylinders labelled $i=1,\ldots,N$, with orientations defined by unit vectors ${\bf u}_{i}$ with components $u_{i\alpha}$, the traceless symmetric nematic order parameter $\bf Q$ has components $Q_{\alpha\beta}=({1}/{N})\sum_{i=1}^{N}\left[({3}/{2})u_{i\alpha}u_{i\beta}-\dfrac{1}{2}\delta_{\alpha\beta}\right]$. The scalar nematic order parameter $S$ is the largest eigenvalue of $\bf Q$.

Hereafter we work in reduced units defined in terms of the system parameters $\epsilon$ and $D$: temperature $T^{*}=k_{B}T/\epsilon$, pressure $P^{*}=Pv_{hsc}/(k_{B}T)$, packing fraction $\eta=v_{hsc}\rho$, where $\rho=N/V$ and $v_{hsc}=\pi D^{2}(D/6+L/4)$ is the volume of a spherocylinder.

Activity in our system is introduced by connecting half of the particles to a thermostat of higher temperature, while maintaining the temperature of the other half fixed at a lower value. Let $T_{h}^{*}$ and $T_{c}^{*}$ be the temperatures of the baths connected to the hot and cold particles respectively, controlled by a Berendsen thermostat Berendsen et al. (1984) with a time constant $\tau_{T}=0.01$. We then define the activity $\chi=(T_{h}^{*}-T_{c}^{*})/T_{c}^{*}$. Starting from a statistically isotropic structure at a definite temperature with $T_{h}^{*}=T_{c}^{*}=5$, we gradually increase the temperature of the hot particles $T_{h}^{*}$, keeping the volume of the simulation box constant (NVT ensemble) Chattopadhyay et al. (2021).

In equilibrium, we observe four stable phases for $L/D=5$: crystal (K), smectic-A (SmA), nematic (N), isotropic (I); three stable phases for $L/D=3$: crystal, smectic-A, isotropic and two stable phases for $L/D=2$: crystal, isotropic (Fig. S1 in Supplemental Information (SI) sp ). Our results are consistent with those of a previous study carried out by Cuetos et al. Cuetos et al. (2002); Cuetos and Martínez-Haya (2015) for soft rods.

To define a criterion analogous to that of Onsager Onsager (1949) for our case, we construct an effective hard-cylinder diameter for the SRS, in terms of the interaction potential and the temperature of the system [Fig. S2], as in Cuetos et al. Cuetos et al. (2005); Cuetos and Martínez-Haya (2015); Boublík (1976):

$$D_{eff}(T)=\int_{0}^{\infty}(1-\exp[-\beta U_{SRS}(d_{m})]dd_{m}.$$

(2)

Therefore, a SRS with aspect ratio $A_{SRS}=L/D$ can be mapped to a HSC with an effective aspect ratio $A_{HSC}=L/D_{eff}$. In Table 1, we mention different values of $A_{SRS}$ and the corresponding values of $A_{HSC}$ at $T^{*}=5$. Using Eq. 2, the value of $A_{SRS}$ corresponding to $A_{HSC}=3.70$ (Onsager’s limit for HSC) becomes $A_{SRS}=3.52$ (Onsager’s limit for SRS) at $T^{*}=5$. This approximate version of Onsager’s criterion is verified by the absence of a nematic phase for SRS at thermal equilibrium with $A_{SRS}=3$ and $2$ in our simulations and those of Cuetos et al. Cuetos et al. (2002); Cuetos and Martínez-Haya (2015).

Starting from a homogeneous isotropic structure, we observe local phase separation between hot and cold particles, which emerges at the macroscopic scale by forming a well-defined interface (Fig. 1). The phase separation is quantified through the difference between the local densities of hot and cold particles. Dividing the simulation box into sub-cells labelled $i=1,...,N_{cell}$ and letting $n^{i}_{hot}$, $n^{i}_{cold}$ be the numbers of hot and cold particles respectively in each cell, we define $\phi=N_{cell}^{-1}\langle\sum_{i=1}^{N_{cell}}{\lvert(n^{i}_{hot}-n^{i}_{cold})\rvert}/{n^{i}_{tot}}\rangle$ Chari et al. (2019); Chattopadhyay et al. (2021) where the average $\langle...\rangle$ is carried out over a sufficiently large number of steady-state configurations. The choice of $N_{cell}$ depends on the aspect ratio of the rod. It is chosen such that each cell contains a sufficient number of particles to get stable statistics.

In Fig. 2-(a), we plot $\phi$ as a function of activity $\chi$ for $L/D=5,3$ and $2$. For the sake of comparison, the system is chosen with a packing fraction between $\eta=0.33-0.36$ for which the system is in the isotropic phase at thermal equilibrium for the given $L/D$ at $T^{*}=5$ (Fig. S1). From Fig. 2-(a), we observe: (i) phase separation starts at a lower activity for higher aspect ratios; (ii) the amount of phase separation at a given $\chi$ is higher for higher aspect ratios. To understand observation (i) precisely, we calculate the critical activity $\chi_{c}$, which is defined as the value of $\chi$ above which macroscopic phase separation is seen (see SI and Fig. S3 for details sp ). The calculated ranges of critical activities for the given packing fractions are: $\chi_{c}=1.4-2$ for $L/D=5$, $\chi_{c}=2-3$ for $L/D=3$ and $\chi_{c}=3-5$ for $L/D=2$.

After phase separation, the interactions between hot and cold particles mostly take place at the interface. The cold zone undergoes an ordering transition above a critical activity $\chi^{*}$ that depends on the aspect ratio of the rods as well as their packing fractions. In Fig. 2-(b), we see that the nematic order parameter of cold particles $S_{cold}$ increases with $\chi$ for $L/D=5$ and $3$ for $0\leq\chi\leq 9$, while for $L/D=2$ it increases above $\chi\geq 10$ (Fig. S4) resulting in a higher $\chi^{*}$ . However, the difference between the critical activities for phase separation and ordering decreases with increasing density. In Fig 2-(c), we have plotted $S_{cold}$ versus $\chi$ for $L/D=2$ at a slightly higher packing fraction, $\eta=0.45$, which also corresponds to the isotropic phase in equilibrium. Here we see both phase separation and ordering transition for the given range of activities. The critical activities for phase separation and ordering at $\eta=0.45$ are $\chi_{c}=3-5$ and $\chi^{*}=5$, respectively.

The phases of the ordered structures in the cold zones are identified by calculating the local nematic order parameter $S$ and appropriate pair correlation functions (see SI for details). We observe that with increasing activity both translational and orientational correlations in the cold zone are enhanced, indicating incipient order in the cold zone. Interestingly, we observe liquid crystal phases for small values of $L/D$ that do not occur for the same parameter range in a one-temperature system, i.e., in equilibrium.

In Fig. 3 (a-d), we plot different pair correlation functions of the cold particles at different activities for $L/D=3$ at $\eta=0.33$. We see that, both of translational $[g(r)]$ and orientational $[g_{2}(r)]$ pair correlation functions are flat in absence of activity ($\chi=0$) which is obvious for an isotropic phase. $g(r)$ develops the first peak at $\chi=4$ and eventually the other peaks at higher values of $\chi$. In Fig. 3(b), we see that, $g_{2}(r)$ has a finite correlation length of roughly $2.5D$, beyond which it decays to zero for $\chi=4,5$ and to a finite value for $\chi\geq 6$. This is also observe in the calculation of the half width at half maximum (HWHM) of $g_{2}(r)$ defined as the distance from the first peak at which the value of $g_{2}(r)$ is half of its value at the first peak [see Fig. S5 in SI for details]. These observations suggest that there exists a finite orientational order in the cold zone for $\chi=4,5$ which designates this phase as nematic. Above this activity, $g_{2}(r)$ develops multiple peaks at longer distances and saturates at a finite value. This is due to the presence of multiple clusters of different average directors, which effectively suppress the overall orientational correlation. In these cases, we calculate the $g_{2}(r)$ in a single cluster of a definite director and find it to saturate at a higher value, as shown in the inset of Fig. 3(b). Smectic and crystalline structures are identified by calculating translational correlations along the parallel [$g_{\parallel}(r)$] and perpendicular [$g_{\perp}(r)$] directions of the average nematic director of the spherocylinders. The periodic oscillations in $g_{\parallel}(r)$ and liquid-like structure in $g_{\perp}(r)$ at $\chi=6,7$ indicate that the phase is smectic, as shown in Fig. 3 (c, d).

Similarly, we find that the system with $L/D=2$ exhibits smectic ordering at $\eta=0.45$, $\chi=9$ as shown in Fig. S6. However, the hot zone shows isotropic structure with reduced packing fractions for each of the cases. The snapshots of the different liquid crystal phases for the respective $L/D$s are shown in Fig. 1. We summarize these results in Table 1 where we mention the possible phases in equilibrium and the emergent phases in active systems for the respective aspect ratios.

To find the microscopic origin of the ordering transitions that are not observed in equilibrium, we calculate the local pressure in the phase-separated system. We divide the simulation box into a number of slabs ($i$) along the direction normal to the interface, following the procedure mentioned in Ref. Chattopadhyay et al. (2021). The location of the interface is identified by calculating the local density of the hot and cold particles. The region where the local densities change sharply between their values in the segregated zones is identified as the interfacial region [Fig. 4(a)].

We then calculate the pressure components along the normal and tangential directions of the interfacial plane using the diagonal components of the stress tensor. We designate the normal direction of the interface as $x$ and the two tangential directions of the interfacial plane as $y$ and $z$. Therefore, the normal $P_{n}$ and tangential $P_{t}$ components of the pressure are defined as: $P_{n}(i)=P_{xx}(i)$ and $P_{t}(i)=(P_{yy}(i)+P_{zz}(i))/2$, where,

$$P_{\alpha\beta}(i)=\frac{1}{V(i)}\left(\sum_{j=1}^{n(i)}mv_{j}^{\alpha}v_{j}^{\beta}+\sum_{j=1}^{n(i)-1}\sum_{k>j}r_{jk}^{\alpha}f_{jk}^{\beta}\right).$$

(3)

The $1^{st}$ and $2^{nd}$ term in Eq. 3 represent the kinetic and virial contributions (arises due to the particle’s interaction) of the pressure tensor, where, $\vec{r_{j}}=(r_{j}^{\alpha})$ and $\vec{v_{j}}=(v_{j}^{\alpha})$ are the position and velocity of the $j^{th}$ particle $(\alpha=x,y,z)$ and $\vec{r_{jk}}$ and $\vec{f_{jk}}$ are the relative distance and interacting force between the SRSs $j$ and $k$. $n(i)$ and $V(i)$ represent total number of particles and volume of the $i^{th}$ slab.

In Figs. 4(c), for the aspect ratio $L/D=3$, we find that, (i) in equilibrium ($\chi=0$), normal and tangential pressure components are equal, implying that the pressure tensor is isotropic throughout the simulation box. (ii) After phase separation ($\chi\geq\chi_{c}$), $P_{n}$ and $P_{t}$ are equal (within the error bars) in the hot zone. Thus, the pressure anisotropy is roughly equal to zero in this region. (iii) At the interface, the tangential pressure decreases and acquire a lower value in the cold zone. However, the normal pressure remains balanced throughout the simulation box. This causes a pressure anisotropy at the hot-cold interface that persists in the bulk cold zone as well [Fig.4(d)]. The pressure anisotropy increases with activity (see Fig. S7) and acquires a finite positive value in the cold zone. This causes an effective compression of the cold zone along the normal direction of the interface that may drive the cold particles to orient along the parallel direction of the interfacial plane, thereby inducing an ordering transition.

To understand further the effects of the activity, we calculate the local heat flux $\vec{J}$ in each slab using the following equation Hansen and McDonald (2013); Bhuyan et al. (2020)

$$\vec{J}(i)=\dfrac{1}{V(i)}\langle\sum_{j=1}^{n(i)}\vec{v_{j}}e_{j}+\sum_{j=1}^{n(i)}\mathbf{\sigma_{j}}\cdot\vec{v_{j}}\rangle,$$

(4)

where $e_{j}=({1}/{2})mv_{j}^{2}+\sum_{j\neq k}U_{jk}$ is the total energy and $\mathbf{\sigma}_{j}=({1}/{2})\sum_{j\neq k}\vec{r_{jk}}\vec{f_{jk}}$ is the stress tensor. In Fig. 4(b), we show the spatial variation of $\vec{J}$ along the direction normal to the interface. We find $\vec{J}=0$ in equilibrium, but it obtains a finite value in the phase separated systems, both at the interface and in the bulk region. This reveals that, though the interaction between hot and cold particles takes place mainly at the interface, its effect extends to the bulk region as well. The sign of $\vec{J}$ indicates heat flows from the bulk hot to the bulk cold zone. This results in heterogeneous activity and broken time-reversal invariance throughout the simulation box giving rise to anomalous thermodynamic behavior (such as pressure anisotropy) away from the interface.

To check for possible finite size effects, we have done a similar analysis for larger system sizes with $N=4608$ for $L/D=2$ and $N=3456$ for $L/D=3$. We observe a similar trend of ordering transition for the respective aspect ratios. Here also we observe a nematic phase for $L/D=3$ and a smectic phase for $L/D=2$ (See SI for details, Fig. S8-S11). Further increase of activity causes the cold zone to crystallize.

Our simulation study examines the effect of two-temperature activity in a soft-rod fluid for a range of effective aspect ratios of the rods. We show that the two-temperature model can give rise to liquid crystal phases that are not observed in equilibrium for the respective aspect ratios. We observe a smectic phase for $L/D=2$ and a nematic phase for $L/D=3$. We find that the presence of two temperatures causes a pressure anisotropy extending from the hot-cold interface into the bulk of the cold zone, and a heat current flowing from the hot to the cold zone. Thus, the nonequilibrium behavior is not limited to the hot-cold interfaces but pervades the system as a whole, driving the anomalous ordering transitions in the cold zone. An understanding of these results within analytical theory, experimental realizations of two-temperature systems, presumably in suspensions with bidisperse motility, and methods to capture and stabilize the anomalously ordered domains are some of the challenges that emerge from our work.

We would like to thank Prof. Aparna Baskaran for helpful discussions. We thank SERB, India for financial support through providing computational facility. JC acknowledges support through an INSPIRE fellowship. JC thanks Tarun Maity, Subhadeep Dasgupta and Swapnil Kole for insightful discussions. SR was supported by a J C Bose Fellowship of the SERB, India, and by the Tata Education and Development Trust, and acknowledges discussions during the KITP 2020 online program on Symmetry, Thermodynamics and Topology in Active Matter. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. CD was supported by a Distinguished Fellowship of the SERB, India.