Orientational ordering of closely packed Janus particles

TL;DR

This study uses Monte Carlo simulations to explore how Janus particles in 2D pack and order, revealing a phase transition to stripe patterns with universality similar to the 3-state Potts model, and analyzing the dynamics of domain growth.

Contribution

It demonstrates a phase transition in 2D Janus particles from disordered to stripe-ordered phases and identifies its universality class, providing insights into their orientational ordering.

Findings

Disordered plastic crystal phase at small patch sizes.

Stripe phase with broken three-fold symmetry at larger patches.

Second-order phase transition with 3-state Potts universality.

Abstract

We study orientational ordering of $2$-dimensional closely packed Janus particles by extensive Monte Carlo simulations. For smaller patch sizes the system remains in the plastic crystal phase where the rotational degrees of freedom are disordered down to the lowest temperatures. There the liquid consist of dimers and trimers of the attractive patches. For large enough patch sizes, the system exhibits a thermodynamic transition into a phase with stripe patterns of the patches breaking the three-fold rotational symmetry. Our results strongly suggests that the latter is a 2nd order phase transition whose universality is the same as that of the $3$-state Potts model in $2$-dimensions. Furthermore we analyzed the relaxation dynamics of the system performing quenching simulations into the stripe phase. We found growing domains of the stripes. The relaxation of key dynamical quantities follow universal scaling features in terms of the domain size.

Figures

Equations

$$H=\sum _{<i,j>}V({\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{S}}_i,{\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{S}}_j)$$

$$V(x,y)=-ϵ\Theta (x-\cos {\theta }_0)\Theta (-y-\cos {\theta }_0)$$

$$C/N{k}_{\mathrm{B}}=({\beta }^2/N)(⟨E^2{⟩}_{\mathrm{eq}}-⟨E{⟩}_{\mathrm{eq}}^2)$$

$$S(\mathbf{k})=\frac{1}{N}\left|\sum _{i=1}^N{e}^{i{\mathbf{r}}_i\cdot \mathbf{k}}\mathbf{S}(0)\cdot \mathbf{S}({\mathbf{r}}_i)\right|,$$

$$\mathbf{k}=\frac{n_1}{L}{\mathbf{b}}_1+\frac{n_2}{L}{\mathbf{b}}_2(n_1,n_2=0,1,2,...,L-1),$$

$${\psi }_x=\frac{\sqrt{3}}{2}(S({\mathbf{k}}_2)-S({\mathbf{k}}_3)),$$

$${\psi }_y=S({\mathbf{k}}_1)-\frac{1}{2}(S({\mathbf{k}}_2)+S({\mathbf{k}}_3)).$$

$$B=\frac{3}{2}\left(1-\frac{⟨{\psi }^4{⟩}_{\mathrm{eq}}}{3⟨{\psi }^2{⟩}_{\mathrm{eq}}^2}\right).$$

$$\varphi (t)=1-\frac{1}{N}\sum _{i=1}^N{\delta }_{{\sigma }_i(t),0}.$$

$$C_{\mathrm{rh}}(r,t)=P(r,t)-\underset{r^{\prime }\to \infty }{\lim }P(r^{\prime },t)$$

$$P(r,t)=\frac{{\sum }_{i,j=1}^N{\delta }_{r,r_{ij}}{P}_{ij}(t)}{{\sum }_{i,j=1}^N{\delta }_{r,r_{ij}}}{P}_{ij}(t)={\delta }_{{\sigma }_i(t),{\sigma }_j(t)}(1-{\delta }_{{\sigma }_i(t),0}{\delta }_{{\sigma }_i(t),0}),$$

$$C_{\mathrm{rh}}(r,t)=\tilde{C}\left(\frac{r}{R(t)}\right).$$

$$\Delta E=E-E_{\mathrm{GS}}\sim R(t{)}^{d_s-d}{d}_{\mathrm{s}}-d=-1$$

$$1-\varphi (t)\sim R(t{)}^{d_s-d}{d}_{\mathrm{s}}-d=-1,$$

$$C(t,t_{\mathrm{w}})=\frac{1}{N}\sum _{i=1}^N{\mathbf{S}}_i(t)\cdot {\mathbf{S}}_i(t_{\mathrm{w}}).$$

$$C(t,t_{\mathrm{w}})=f\left(\frac{R(t)}{R(t_{\mathrm{w}})}\right)\sim {\left(\frac{R(t)}{R(t_{\mathrm{w}})}\right)}^{-\lambda }(t\gg t_{\mathrm{w}}).$$

Full text

I Introduction

Colloids offer possibilities to create new condensed matters starting from the design of the building blocks: the size, shape and mutual interactions of the colloidal particles. The primary interaction between the colloids is the hard-core repulsion. Oftenly they also have attractive interactions on top of it. If the attractive interaction is sufficeintly short-ranged, they are called as sticky colloids which attract each other if and only if they are in contact: their surfaces are separated only by a negligibly short distance. Sticky colloids exhibit various interesting phenomena such as crystal-crystal transitionsBolhuis and Frenkel (1994); Shin and Schweizer (2013), glass-glass transitionsFabbian et al. (1999); Sciortino (2002), gas-liquid critical phenomenaMiller and Frenkel (2003), gelationLu et al. (2008), etc.

Moreover colloids can have rotational degree of freedoms. If the colloids are non-spherical like ellipsoids, the mutual interaction becomes dependent not only on the mutual distance between the centers of the particles but also on their orientations in complex manners. The resultant anisotropic interactions bring about richer variation of gels Madivala et al. (2009); Ruzicka et al. (2011). The so called patchy colloids, which have the sticky attractive interactions but only through their patch like limited regions on their surfaces, provide varieties of soft condensed matters. In particular simple spherical colloids with designed patches on top of them offer opportunities to study various non-trivial consequences of anisotropic interactions in a controlled way. It is known that patchy colloids realize various stable phases depending on the interaction strength, number and size of the patches.Bianchi et al. (2006); Sciortino et al. (2009, 2010); Romano and Sciortino (2011); Chen et al. (2011); Iwashita and Kimura (2013); Vissers et al. (2013)

In the present paper, we study the simplest spherical patchy colloids with just one patch on each of them, which are called as ”Janus particles”, put in the $2$-dimensional closest packing configuration, i.e. the triangular lattice of the spherical particles directly touching all the nearest neighbors. (see Fig. 1(a)). In such an configuration particles cannot move but just rotate. Thus the configuration of the system can be described solely by ’spins’ or ’directors’ $\bm{S}_{i}$ ($i=1,2,\ldots,N$) which represent the orientations of the patches of the particles $i=1,2,\dots,N$. This is a very interesting set up which realizes a simplest plastic crystal: orientationally disordered states on top of a translationally ordered, crystalline background. Indeed such set ups have been implemented and analyzed experimentally.Iwashita and Kimura (2014); Jiang et al. (2014). However the nature of the orientational ordering have not been fully studied yet. While a theoretical study Shin and Schweizer (2014) predicts an orientationally ordered state for large enough patch sizes, an experiment Jiang et al. (2014) suggests an orientational glass phase. In the present work, we perform extensive Monte Carlo simulation to study both the static and dynamic aspects of the orientational configurations.

The attractive interaction between two Janus particles is switched on only if their patches are directly touching with each other as in Fig. 1 (a). Let us say that the two particles have an attractive contact in such a situation. Then a key parameter of the system is the patch size which can be parameterized by the angle $\theta_{0}$ (see Fig. 1 (b)). The spatial patterns of the directors in the ground states changes dramatically with the patch size (see Fig. 2 (a)).

1. Dimers ($0^{\circ}<\theta_{0}\leq 30^{\circ}$): Each particle can have attractive contacts at most with 1 other adjacent particle so that a ground state configuration consists of dimers of patches. This is just a dimer covering problem on the triangular lattice so that there are energetically degenerate, numerous ground states.

2. Trimers ($30^{\circ}<\theta_{0}\leq 60^{\circ}$) Each particle can have attractive contacts at most with 2 other particles so that a ground state configuration consists of trimers of patches. Again energetically degenerate, numerous disordered ground states can be expected. 3) Stripes ($60^{\circ}<\theta_{0}\leq 90^{\circ}$): Each particle can have attractive contacts at most with 3 other adjacent particles. A ground state consists of stripes of the patches.Shin and Schweizer (2014)

One would naturally expect rotationally disordered states (plastic crystals) with dimers or trimers at low enough temperatures in the regime 1) and 2). A possibility in the regime 3) is a long-ranged ordering of the stripes. Shin and Schweizer studied and analyzed the system theoretically based on a harmonic approximation from the ground states and proposed a theoretical phase diagram. Shin and Schweizer (2014) They argued that there are dimer liquid, trimer liquid and crystalline stripe phases at low temperatures which are separated from the high temperature, totally disordered orientational liquid (plastic crystal) phase by 1st order transitions. However the validity of the harmonic approximation up to higher temperatures may be questioned. Thus it is not obvious whether all the transitions are real thermodynamic transitions nor they are all first order transitions. Rotational ordering of the closely packed Janus particles have also attracted interests in experimentsIwashita and Kimura (2014) which found qualitative agreements with the theoretical phase diagram but the precise nature of the transitions remains unknown. In order to clarify these issues we performed a detailed analysis of the putative orientational phase transitions by equilibrium Monte Carlo simulations.

Dynamics of the rotational degree of freedom was also studied experimentally in Jiang et al. (2014) which reported that Janus particles with patches size $\theta_{0}=90^{\circ}$, making the stripe patterns, show glassy dynamics avoiding long-ranged ordering of stripes. The possibility of creating an orientational glass from a plastic crystal is a very interesting question by itself Suga and Seki (1974) which deserves to be investigated thoroughly. Indeed if the equilibrium transition to the stripe phase is a 1st order transition as suggested theoretically, super-cooling of the orientationally disordered state and the ultimate orientational glass transition is not inconceivable. In order to clarify this issues we performed a detailed analysis of the orientational dynamics after temperature quenches by Monte Carlo simulations.

In this paper, we report the following two sets of main results. The first is from the equilibrium Monte Carlo simulations. We found that the transition to the dimer and trimer states are not thermodynamic phase transitions but smooth crossovers. Moreover we found that transition to the stripe phase is not a 1st order transition but a 2nd order phase transition. In the stripe phase there are 3 possible directions of the stripes and the phase transition breaks the $3$-fold rotational symmetry. We defined an appropriate order parameter which is able to detect the breaking of the 3-fold rotational symmetry. We found that the universality class of the transition is the same as $3$-state Potts model in 2 dimensions. The second set of results is obtained from the analysis on the relaxational dynamics of the system after temperature quenches by the Monte Carlo simulations. We found that the stripe phase ordering kinetics obey a universal scaling feature of the standard coarsening phenomena instead of the glassy dynamics at odds with the experimental suggestion.

The organization of the paper as follows. In sec. II we introduce the model and in sec. III we explain our simulation method. Then in sec. IV we report the results of equilibrium Monte Carlo simulations. And in sec. V we report the results of out-of equilibrium Monte Carlo simulations. In sec. VI we discuss implications of our results for experiments. Finally in sec .VII we conclude the paper with a few remarks.

II Model

In the present paper we focus on the rotational degree of freedoms of the closely packed Janus particles in 2 dimension (Fig. 1 (a)) freezing out the translational motions for simplicity. The center of mass position ${\bf r}_{i}$ of the Janus particles $i=1,2,\ldots,N$ are put on the triangular lattice whose three ridges being parallel to the following three unit vectors ${\bf a}_{1}=(1,0)$, ${\bf a}_{2}=(1/2,\sqrt{3}/2)$, ${\bf a}_{3}=(-1/2,\sqrt{3}/2)$. Then the system can be considered as an assembly of classical vectors representing the directors $\bm{S}_{i}$ ($i=1,2,...N$) put on the triangular lattice. For simplicity we describe the directors by $2$-component spin vectors $\bm{S}_{i}=(S^{x}_{i},S^{y}_{i})$ which lie in the same plane of the triangular lattice. Note that here we are disregarding the out-of-plane components of the directors which we believe to be irrelevant. The Hamiltonian is given by,

[TABLE]

with

[TABLE]

which is the rotational part from Kern-Frenkel potential. Kern and Frenkel (2003) Here $\hat{\bm{r}}_{ij}=({\bf r}_{j}-{\bf r}_{i})/|{\bf r}_{j}-{\bf r}_{i}|$ and $\Theta(x)$ is the Heaviside’s step function : $\Theta(x)=1(x\geq 0),-1(x<0)$. Note that on the triangular lattice, $\hat{\bm{r}}_{ij}$ is parallel to one of the three vectors ${\bf a}_{1}$,${\bf a}_{2}$,${\bf a}_{3}$ shown in Fig. 1 (a). Patches size is parameterized by the angle $\theta_{0}$.

It can be easily seen that the potential $V$ takes a constant negative value $-\epsilon$ ($\epsilon>0$) if and only if the two particles are in an attractive contact, i.e. they are in contact through their patches. As the result, the energy levels of the system become discretized in spite of the fact that the spins are continuous variables. We introduce a dimensionless energy $\tilde{E}=H/\epsilon$ (hereafter we call it simply as $E$) and a dimensionless temperature $\tilde{T}=k_{B}T/\epsilon$ (hereafter we call it simply as $T$) where $k_{\rm B}$ is the Boltzmann’s constant. We denote the thermal average in equilibrium as $\langle\ldots\rangle_{\rm eq}$.

Let us comment on the symmetry of the system. The system is invariant under global rotations, by which not only the spins but also the whole lattice are rotated all together. Note that the system is not invariant under rotations of the spins only at variance with usual spin systems. Since we have fixed the translational degree of freedoms in the present study, the global rotational symmetry is anyway lost. More important is the fact that the spin Hamiltonian is invariant global rotations by a discrete angle $2\pi/3$ reflecting the triangular lattice. The ground states with the ordered stripe patterns (see Fig. 2 (a)) break the 3-fold rotational symmetry.

III Monte Carlo simulation

We performed Monte Carlo simulations with the Metropolis method.Landau and Binder (2014) In order to facilitate the computation using vector computers, we employed the sub-lattice flipping for spin updates. To update a spin $\bm{S}_{i}$, a new candidate spin configuration $\bm{S}^{\prime}_{i}$ was created by rotating $\bm{S}_{i}$ by a random angle $\gamma$, which obeys a uniform distribution in the range $-\pi/3\leq\gamma<\pi/3$. In a unit of Monte Carlo step (MCS) a spin update is tried for all spins $\bm{S}_{i}$ ($i=1,2,\ldots,N$). In the present work, we used system size $L=48,72,108$ for the equilibrium simulations and $L=144,288,576,1152$ for the out-of equilibrium simulations. We employed periodic boundary conditions in both simulations.

• •

Equilibrium Monte Carlo simulations

In equilibrium simulations, we observed static physical quantities for various patches sizes $\theta_{0}$ and different temperature $T$. The equilibrium states are realized by following slow cooling procedure. (0) Fully equilibrated configurations are prepared at a high enough temperature $T=0.70$ for each patches size $\theta_{0}$. (1) The temperature is lowered by a small step $\Delta T=0.01$. (2) Then the system is thermalized by $3\times 10^{5}$ (MCS). (3) Observations of physical quantities are made during additional $3\times 10^{5}$ (MCS) to evaluate their thermal averages $\langle\ldots\rangle_{\rm eq}$ by their time averages. (4) return to (1).

In the vicinity of the transition temperatures, the temperature step $\Delta T=0.01$ used in the annealing is replaced by a smaller one $\Delta T=0.002$. To conform that the equilibration time is large enough, we performed simulations using a longer equilibrating time $3\times 10^{6}$ (MCS) and checked that the results do not change.

We performed $50$ statistically independent runs and evaluated the time averages and the mean-squared errors of the physical quantities.

• •

Out-of equilibrium Monte Carlo simulations

To study the relaxation dynamics after temperature quenches into the stripe phase, we mainly worked with the patches size $\theta_{0}=90^{\circ}$. The unit of time is $1$ (MCS). We quenched the system from $T=\infty$ to $T=0$ at $t=0$ (MCS). More precisely, we took totally random spin configurations $\bm{S}_{i}$ as the initial spin configuration and performed Monte Carlo simulations at $T=0$: each spin flip is accepted only if it lowers the energy of the system. To check possible dependence on the details of the algorithm used for the spin updates, we compared results with ones obtained by single flipping instead of the sub-lattice flipping.

We performed $120$ statistically independent runs for $L=144,288$, $90$ ones for $L=576$, $60$ ones for $L=1152$ and took averages of all the physical quantities over the runs and also estimated their mean-squared errors.

IV Equilibrium Monte Carlo simulations

IV.1 Phase diagram

First let us examine the phase diagram. Closely packed Janus particles exhibit three different types of ground states: dimer $0^{\circ}<\theta_{0}\leq 30^{\circ}$, trimer $30^{\circ}<\theta_{0}\leq 60^{\circ}$ and stripe states $60^{\circ}<\theta_{0}\leq 90^{\circ}$ Shin and Schweizer (2014); Iwashita and Kimura (2014). It has been proposed theoretically that there are three low temperature phases corresponding to the three distinct types of ground states. In the theory Shin and Schweizer (2014) the dimer and trimer phases are considered as liquid phases while the stripe phase is considered as a crystalline phase. Thus the 1st basic issue to clarify is actual existence of such thermodynamic phases at low temperatures. We must note however that in the present paper we are not taking into account the vibrational motions of the particles which is included in the theory Shin and Schweizer (2014). We assume that the vibrational motions do not significantly change the nature of the phases.

In Fig. 2 (a) we show the phase diagram we determined based on the measurement of the heat-capacity as we explain below. There we also show the spatial patterns of the ground states mentioned above. We define the heat capacity $C$ as,

[TABLE]

where we used the inverse temperature $\beta=1/k_{\rm B}T$. Here after we call $C/Nk_{\rm B}$ simply as $C$. We measured it by performing equilibrium Monte Carlo simulations at various patch sizes $\theta_{0}$ and temperatures $T$ (more precisely $\tilde{T}=k_{\rm B}T/\epsilon$). At each values of $\theta_{0}$ we find a peak of the heat capacity at a temperature $T_{\rm peak}(\theta_{0})$ as shown in Fig.2 (b). It can be seen that $T_{\rm peak}(\theta_{0})$ slightly depend on the system size $L$ but only very weakly. In the panel (a) we display the line $1/T_{\rm peak}(\theta_{0})$ vs $\theta_{0}$ obtained in the system of size $L=72$. The qualitative feature of the peak temperature $T_{\rm peak}(\theta_{0})$ appears indeed similar to the phase boundary between the high temperature plastic phase and low temperature phases proposed theoretically Shin and Schweizer (2014) (see Fig. 7 (a),(b) of the latter reference), which was observed also by a numerical analysis of the spatial patterns of the configurations Iwashita and Kimura (2014) (See Fig. 10 (a) of the latter reference).

In the range $0^{\circ}<\theta_{0}\leq 60^{\circ}$ we found no singularity at $T_{\rm peak}(\theta_{0})$ but just a broad peak as shown in Fig. 2 b). Furthermore we found are no appreciable system size dependencies there. These results suggest $T_{\rm peak}(\theta_{0})$ should be regarded just as a crossover line and that the dimer and trimer phases are not true thermodynamics phases. On the contrarily, in the range $\theta_{0}>60^{\circ}$ we find a singular behavior of the heat-capacity at $T_{\rm peak}(\theta_{0})$ as shown in panel b) suggesting a thermodynamic phase transition.

IV.2 Stripe order

We have seen that the heat-capacity exhibits a singular behavior suggesting a thermodynamic phase transition for $\theta_{0}>60^{\circ}$ into the stripe phase. Then an important question is whether it is a 1st order transition as suggested theoretically Shin and Schweizer (2014) or a continuous transition. In order to examine the nature of the low temperature phase and the phase transition we performed further analysis as the following. In the following we only display results for $\theta_{0}=90^{\circ}$ for clarity but similar results were obtained at other values of $\theta_{0}(>60^{\circ})$ as well.

In order to characterize the low temperature phase, we examined the structure factor defined as,

[TABLE]

where $\bm{k}$ is a wave vector and $\bm{b}_{1}=2\pi(1,1/\sqrt{3}),\bm{b}_{2}=2\pi(0,2/\sqrt{3})$ are reciprocal vectors of the triangular lattice. We display in the inset of Fig.  2 c) an intensity plot of the structure factor which reveals prominent three distinct spots at $\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}$ ($(n_{1}/L,n_{2}/L)=(1/2,0),(0,1/2),(1/2,1/2)$ indicating the three possible orientations of the stripes.

If a stripe structure exhibits a long-ranged order in the low temperature phase, it means breaking of the 3-fold rotational symmetry. In order to detect the possible spontaneous symmetry breaking, let us introduce an order parameter $\bm{\psi}=(\psi_{x},\psi_{y})$ defined by,

[TABLE]

Clearly, it will be averaged out to [math] in a thermal average if the $3$-fold rotational symmetry is preserved. In Fig. 2 (c) we show the temperature dependence of the absolute value of the order parameter $\langle|\bm{\psi}|\rangle_{\rm eq}$. It can be seen that it rapidly increases in the vicinity of $T_{\rm peak}(\theta_{0}=90^{\circ})$ by lowering the temperature strongly suggesting spontaneous breaking of the 3-fold rotational symmetry.

To analyze the phase transition more precisely, we examined the Binder’s parameter defined by,

[TABLE]

As shown in Fig. 2 (e), the data of different system sizes exhibit a crossing feature strongly suggesting a continuous phase transition. From the latter we estimate the critical temperature of the system as $T_{c}\approx 0.275$ at $\theta_{0}=90^{\circ}$, which is in good agreement with the peak temperature $T_{\rm peak}(\theta_{0}=90^{\circ})$ of the heat-capacity (See Fig. 2).

Finally, to examine the universality class of the continuous phase transition, we performed finite size scaling analyses as follows. Since the order parameter has the $3$-fold rotational symmetry, it is natural to expect that the universality class of the transition is the same as that of the 3-state Potts model in 2d, whose critical exponents are known to be $\nu=5/6,\beta=1/9,\gamma=13/9$ Wu (1982). In Fig. 2 (d) and (f) we display finite size scaling plots of the order parameter $\langle|\bm{\psi}|\rangle_{\rm eq}$, its 2nd moment $\langle\bm{\psi}^{2}\rangle_{\rm eq}$ and the Binder’s parameter $B$. The successful collapse of the data in the scaling plots strongly suggests that the universality class of the phase transition into the stripe phase is indeed that of the 3-state Potts model in 2d. Thus our result is qualitatively different from that of the the harmonic theory Shin and Schweizer (2014) which proposes a 1st order phase transition.

V Out-of equilibrium Monte Carlo simulations

In previous section, we demonstrated that there is a continuous phase transition into the stripe phase for patch size in the range $60^{\circ}<\theta_{0}\leq 90^{\circ}$. Now we turn to study the relaxational dynamics after temperature quenches into the stripe phase. Most natural expectation would be that the relaxation of the system proceeds as a coarsening process Bray (2002) of the stripes. However an experimental study Jiang et al. (2014) has suggested some sort of glassy dynamics.

V.1 Local order parameter

By inspecting the stripe structure of the ground state (see Fig. 2 (a)), we first notice that a stripe consists of a train of tetramers. In a tetramer, four particles are in attractive contacts with each other. A tetramer can be represented by a rhombus whose ridges connect the centers of the particles in the tetramer as shown in Fig. 3 (a). There are three possible orientations of the tetramers (rhombus) due to the 3-fold rotational symmetry of the triangular lattice. We label them by three colors: red, blue and green or the corresponding integers $\sigma=1,2,3$. The latter integers corresponds to vectors $\bf{a}_{\sigma}$ ($\sigma=1,2,3$)(See Fig. 1), each of which is parallel to the vector connecting the shorter diagonal of the corresponding rhombus.

A stripe is a train of tetramers with the same orientation (color) $\sigma$ as shown in Fig. 3 (b). In such a train, two adjacent tetramers are connected by a link parallel to $\bm{a}_{\sigma}$, along which two particles in different tetramers are in an attractive contact. Such a train can be continued indefinitely without any energy cost along two directions among $\bm{a}_{1}$,$\bm{a}_{2}$,$\bm{a}_{3}$ except for $\bm{a}_{\sigma}$.

We performed out-of equilibrium Monte Carlo simulations starting from $T=\infty$ random initial configurations. Since spin updates are performed at $T=0$, the system simply tries to lower the energy increasing the number of attractive contacts. In Fig. 3 (c)-(e) we show some snapshots of the configuration of the system represented by the oriented (colored) rhombuses. There we also display the attractive contacts which are not included in the tetramers as black lines. It can be seen that the length of stripes, namely trains of tetramers with the same orientations, increase with time.

Moreover we find that stripes with the same color also grow into their transverse directions resulting in growth of 2 dimensional domains of tetramers with the same color. This is not so obvious at a first sight because stripes of different colors can run in parallel in a ground state. For instance, stripes of $\sigma=1$ and $\sigma=2$ can be extended in parallel along the same direction $a_{3}$. In other words, there are large number of ground states in which all stripes are directed toward a given direction $a_{\sigma}$ $\sigma\in(1,2,3)$ where each of the stripes can take one of the 2 possible values $\sigma^{\prime}(\neq\sigma)$ arbitrarily. Indeed we can finds some stripes of different colors running in parallel in the snapshots shown in Fig. 3.

Where stripes of different colors running into different directions meet, we find a local region not belonging to tetramers. There the local energy is higher as compared with the bulk of the domains. Thus the latter may be regarded as energetic defects. Energy relaxation should take place by reducing the area of such energetic defects which amount to expand some domains at the expense of the others, i. e. domain growth. Note also that a stripe can be bended without any energy cost as shown in Fig. 3 (b). Indeed we can find such configurations in the bulk of the domains. Presumably this bending mechanism allows formation of 2 dimensional domains made of stripes which are 1 dimensional.

Based on the above observation let us introduce a local order parameter as follows. We define a Potts variable $\sigma_{i}(t)$ at each site $i=1,2,...,N$ such that it equals to the color $\sigma\in(1,2,3)$ of the tetramer (rhombus) it belongs to at time $t$. If it does not belong to any tetramers, we assign $\sigma_{i}(t)=0$.

V.2 Dynamical observable

Now let us introduce several dynamical quantities to study the relaxational dynamics based on the local order parameter $\sigma_{i}(t)$ introduced above. First, we introduce the density $\phi(t)$ of the rhombuses by,

[TABLE]

Note that $1-\phi(t)$ is the density of spins not belonging to the tetramers. We also measure the relaxation of energy density $E(t)/N$.

To characterize the domain growth, a spatial correlation function of the Potts variable can be introduced as,

[TABLE]

with

[TABLE]

where $\delta_{\alpha,\beta}$ is Kronecker’s delta. In simulations, since we can’t take a limit $\lim_{\bm{r}^{\prime}\rightarrow\infty}$, we subtract correlation with the most distant particles, i. e. $r_{ij}=L/2$.

V.3 Domain growth

Now let us discuss our results on the dynamics. In Fig. 4 (a) we display the rhombus spatial correlation function. We can see that the correlation length becomes larger with time reflecting the domain growth. We evaluated the domain size $R(t)$ as the shortest distance $r$ such that $C_{\rm rh}(r,t)<0.01$. As shown Fig. 4 (b), we find that the resultant domain size grows initially as $R(t)\sim t^{0.5}$. However it crossovers to slower growth law $R(t)\sim t^{0.38}$ at larger time scales $t^{\ast}\approx 10^{4}$. The data suggest this crossover is not a finite size effect.

The initial growth law $R(t)\sim\sqrt{t}$ is the same as in the standard coarsening process of non-conserved order parameters Bray (2002). On the other hand, the slower growth law found at larger time scales is not a standard one and deserves to be studied further in the future.

In Fig. 4 (c) we show a scaling plot of the rhombus spatial correlation function assuming a scaling form,

[TABLE]

It is clear that the scaling law is satisfied as in standard coarsening systems Bray (2002). Remarkably the scaling function is universal in spite of the fact that the growth law $R(t)$ itself exhibits the crossover behavior mentioned above. Such robustness of the scaling function of the spatial correlation function is known also in some other cases Bray (2002).

V.4 Energy and rhombus density relaxation

Next let us examine the relaxation of the energy and the rhombus density. Since we are working at $T=0$, the energy is expected to converge to that of the ground state which is $E_{\rm GS}=-1.5N$. A natural scaling ansatz for the relaxation of the energy is

[TABLE]

where $d$ and $d_{\rm s}$ are spatial dimension of the system and fractal dimension of domain wall respectively. The last equation $d_{s}=d-1$ holds if the domain is compact as in usual coarsening systems. In Fig. 4(d) we show the relaxation of the energy density $\Delta E(t)/N$. This result is consistent with the expected scaling law. We checked that the scaling behavior does not change by using single-spin flipping instead of the sub-lattice flipping.

Similarly we expect that the density of the rhombus density behaves as,

[TABLE]

because the density $1-\phi(t)$ which do not belong to the rhombus should correspond to the area where the local energy is high, i. e. domain wall. As shown Fig. 4 (d) this scaling law is also satisfied.

Again, it is remarkable that the scaling behaviors in terms of the domain size $R(t)$ are universal in spite of the fact that the growth law $R(t)$ itself exhibits the crossover behavior as mentioned before.

V.5 Two time autocorrelation function

Finally let us discuss the spin-auto correlation function,

[TABLE]

As shown in Fig. 5 (a), the two time quantity strongly depends on the waiting time $t_{\rm w}$ as expected. A natural scaling ansatz is Fisher and Huse (1988); Bray (2002); Newman and Bray (1990),

[TABLE]

It is clearly satisfied by our data as shown in Fig. 5 (b). The result suggests that the non-equilibrium dynamical exponent $\lambda=1$. Once again, it is remarkable that the scaling behavior in terms of the domain size $R(t)$ is universal in spite of the fact that the growth law $R(t)$ itself exhibits the crossover behavior as mentioned before.

VI Implications for experiments

Let us discuss implications of our results for experiments. In equilibrium simulations discussed in sec.  IV, we made observations by slowly cooling the system (decreasing $k_{B}T/\epsilon$) keeping the system in equilibrium. In experiments of colloidal particles, it may be possible to realize quasi-static, slow changes of the strength of effective interaction $\epsilon$ via changing the properties of the solution Iwashita and Kimura (2013). The structure factor $S(\bm{k})$ is well accessible in experiments Iwashita and Kimura (2014). Thus it would be very interesting to perform systematic analysis of the nature of the possible phase transitions in experiments and make comparisons with our results. We have to note however that the effects of the gravity, friction between the particles and friction between the particle and substrate which may play some roles in experiments are not included in our present model.

The dynamics of the stripes was investigated in a previous work Jiang et al. (2014) both experimentally and numerically. The spin auto-correlation function Eq. (13) was also studied but unfortunately without paying attention to the waiting time $t_{\rm w}$ effect. There the decay of the spin auto-correlation function was fitted to a stretched exponential function $A\exp(-Bt^{\beta})$ which was interpreted as a signature of glassy dynamics. Indeed our data shown in Fig.  5 can also be fitted by such a formula. However the result of our scaling analysis in terms of the domain size $R(t)$ presented in the previous section strongly suggests that the relaxational dynamics, including the aging effect, can be simply understood as a consequence of a coarsening process. Moreover we consider that it is rather unlikely to realize super-cooled, glassy liquid state since it turned out that the transition into the stripe phase is a 2nd order transition rather than a 1st order transition.

Comparing our data with the result presented in Ref Jiang et al. (2014), we expect our 1(MCS) roughly corresponds to $10^{-2}$(sec) in the experiment. Thus we expect the experimental time scale $O(10^{3})$ (sec) is large enough to investigate the domain growth process we reported. However finite size effects should be severe in the previous experiment which uses $L\approx 50$. Indeed we found strong finite size effects for smaller system sizes as shown in Fig. 4 (b).

VII Conclusion

To summarize, we studied the equilibrium and dynamic properties of directors of the closely packed Janus particles by Monte Carlo simulations. Although the gross phase behavior we found is in agreement with the previous theoreticalShin and Schweizer (2014) and experimental works Iwashita and Kimura (2014), the present work clarified the precise nature of the phase transitions. In equilibrium, we found that the transitions from the high temperature plastic phase to low temperature dimer ($0^{\circ}<\theta_{0}\leq 30^{\circ}$) and trimer ($30^{\circ}<\theta_{0}\leq 60^{\circ}$) ’phases’ are smooth crossovers. On the other hand, in $60^{\circ}<\theta_{0}\leq 90^{\circ}$ we found 2nd order phase transition into stripe phase whose universality class turned out to be the same as that of the $3$-state Potts model in $2$-dimensions.

In relaxation dynamics in the stripe phase $60^{\circ}<\theta_{0}\leq 90^{\circ}$ bellow $T_{\rm c}$, we found strong evidences of coarsening process. The growth law of the domain size $R(t)$ exibits a crossover from the usual behavior $R(t)\sim\sqrt{t}$ at shorter time scales to a slower growth at longer time scales. At the moment the mechanism of the crossover is unknown. Nevertheless scaling analysis of various physical observable in terms of the domain size $R(t)$ reveal salient universal features of the coarsening process of stripes. Such robustness of the scaling features in terms of the domain size with some complicated crossover behaviors renormalized into the domain growth law is known in some systems Yoshino et al. (2002); Corberi et al. (2011). We found no evidence of glassy dynamics originally anticipated in Jiang et al. (2014). However we believe the possibility to realize orientational glasses using patchy colloids is a legitimate and interesting problem by itself Yoshino (2017) which should be explored more extensively.

Finally we note that we have not considered at all possible effects of translational motions of the particles in the present paper. In 2 dimensions translational thermal fluctuations can be strong Mermin and Wagner (1966). Also coupling between translational and rotational degree of freedom should bring about more varieties of phenomena. Certainly these points should be examined in the future.

Acknowledgements.

The author thanks Yuliang Jin, Macoto Kikuchi and Takenobu Nakamura for useful discussions. This work was supported by KAKENHI (No. 25103005 “Fluctuation & Structure” and No. 50335337) from MEXT, Japan. This work was achieved through the use of OCTOPUS and supercomputer system SX-ACE at the Cybermedia Center, Osaka University.