Percolation Phase Transition from Ionic Liquids to Ionic Liquid Crystals

Due to their complex molecular structures and interactions, phase behaviors of complex fluids are quite often difficult to be identified by common phase transition analysis methods. Percolation phase transition, on the other hand, only monitors the degree of connection among particles without strict geometric requirements such as translational or orientational order, and thus suitable for pinpointing phase transitions of complex fluids. As typical complex fluids, ionic liquids (ILs) exhibit phases beyond the description of simple liquid theories. In particular, with an intermediate cationic side-chain length, ILs can form the nanoscale segregated liquid (NSL) state, which will eventually transform into the ionic liquid crystal (ILC) structure when the side chains are adequately long. However, the microscopic mechanism of this transformation is still unclear. In this work, by means of coarse-grained molecular dynamics simulation, we show that, with increasing cationic side-chain length, some local pieces of non-polar domains are gradually formed by side chains aligned in parallel inside the NSL phase, before an abrupt percolation phase transition happens when the system transforms into the ILC phase. This work not only identifies that the NSL to ILC phase transition is a critical phenomenon, but also demonstrates the importance of percolation theory to complex fluids.

In contrast to simple liquids, complex fluids are composed of molecules with complex structures and interactions, whose physical properties are beyond the description of simple liquid theories. Particularly, phase behaviors of complex fluids are quite often difficult to be identified by common phase transition analysis methods, such as calculating the system translation or orientation correlation functions. On the other hand, in a condensed matter system, particles often form clusters with various sizes. Before the percolation phase transition, particles aggregate locally to form small separate clusters whose sizes are incomparable to the system size. At the transition point, separate clusters suddenly percolate to form a giant largest cluster whose size is comparable to the system size and the size of the second largest cluster is very small compared to the largest cluster. The percolation phase transition is a critical phenomenon (second-order phase transition) 1 , whose universal critical exponents for various dimensions are listed in Table 1 2 . Percolation phase transition in condensed matter systems has been intensively studied 1,3 , and exhibits its usefulness in glass transition 4-6 and conductivity 7,8 . Nevertheless, percolation phase transition may be an analyzing tool much more powerful beyond the current understanding, especially for the abundant exotic phase behaviors in complex fluids, because it only requires a minimal geometrical feature of connectivity, which is a much looser requirement than spatial symmetries, such as the translational or orientational order.
Ionic liquids (ILs), also known as room-temperature molten salts, are composed of bulky organic cations and smaller anions. Because of their dual ionic and organic nature 9 originated from the subtle balance among various interactions, such as Coloumbic, van der Waals (VDW), hydrogen bonding, etc., ionic liquids are typical complex fluids exhibiting abundant novel phase behaviors simple liquids do not have. In particular, the nanoscale segregated liquid (NSL) state can form when the cationic side chain length is intermediate 10,11 , which transforms into ionic liquid crystals (ILCs) when the cationic side chain is adequately long [12][13][14][15][16][17] . The formation of the NSL state is attributed to the Coulombic interaction between the cationic head groups and anions which forces the cationic side chains to aggregate and form separate nonpolar tail domains. The VDW interaction between cationic side chains increases with side-chain length, and when the side chains are adequately long, the strength of their VDW interaction allows the side chains to align in parallel and form ILCs that have unique features beyond traditional liquid crystals 18,19 , e.g., polar layers are connected by charged atomic groups 20 . Neither phase can form in a simple liquid, so the exact phase transition point between these two phases may not be identified by traditional methods.
In this paper, we perform coarse-grained (CG) molecular dynamics (MD) simulations for [C n MIM][NO 3 ] (n = 6,8, …, 22, abbreviated as C n thereafter) ILs with a simulation size of 4096 ion pairs and an anisotropic barostat allowing the simulated systems to change their sizes in each dimension independently. Our simulation results at T = 400 K indicate that a sharp percolation phase transition appears at n = 18. The corresponding physical picture is that, with increasing cationic side-chain length, the IL structures first gradually change from the globally isotropic NSL state to a liquid state with some local small clusters formed by cationic side chains aligned in parallel, and then go through a sharp percolation phase transition when the majority of the cationic side chains are globally aligned in parallel. Therefore, this IL to ILC phase transition is a critical phenomenon, which is difficult or even impossible to be unambiguously identified by common phase transition analysis methods, demonstrating that percolation phase transition is particularly useful for analyzing phase behaviors of complex fluids.

Results
Structural properties. The [C n MIM][NO 3 ] (n = 12, 14, …, 24) IL systems were modeled with the effective force coarse-graining (EF-CG) force field 21,22 , as shown in Fig. 1. For each system, 4096 ion pairs were put in a parallelepiped simulation box with the periodic boundary condition (PBC) applied to all three dimensions, and the box size in each dimension changes independently in the NPT ensemble. The equilibrium structures at T = 400 K were obtained after appropriate simulated annealing processes (see the Methods section). As shown in Fig. 2a, in agreement with the previous studies 10,11,20 , from the snapshot we can see that C 12 forms the NSL state with a continuous polar network composed of charged groups (anions and cationic head groups) and separate nonpolar domains composed of cationic tail groups. From Fig. 2b, we can see that the C 16 system with 4096 ion pairs has only formed some local pieces of cationic side chains aligned in parallel but without long-range structural order, which cannot be regarded as in the ILC state. By contrast, as can be seen in Fig. 2c, the C 22 system has adequate side chains globally aligned in parallel to form the ILC structure.
To quantify the structural changes, we first calculated the radial distribution functions (RDFs) for the center-of-masses (COMs) of side chains. As shown in Fig. 2d, the first peak of the RDF initially increases slowly from C 12 to C 16 and then increases quickly from C 18 to C 20 , before it reaches a quite high value for C 22 . This result indicates that the aggregation degree of cationic side chains increases almost monotonically with side-chain length, and a drastic change happens at C 18 .
The heterogeneity order parameters (HOPs) 23 were then calculated for anions (CG sites D), cationic head groups (CG sites A), and cationic tail groups (CG sites E) to investigate the change of aggregation degree with side-chain length. As can be seen in Fig. 2e, the HOP values for anions and cationic head groups overall grows gradually with increasing side-chain length, suggesting that the charged groups become more and more  ]. The EF-CG force field was used to model the IL systems at the coarse-grained level, where the cationic imidazolium ring is coarse-grained as CG site A, the single methyl group as CG site B, the tail methyl group as CG site E, CG sites M1, M2, M3 and M4 correspond to the charged methylene groups connected to the imidazolium ring, and CG sites C correspond to all the charge-neutral methylene groups.
aggregated due to increasing VDW interactions among side chains. In contrast to the higher first RDF peak of C 20 than C 18 , the HOP values for the charged groups of C 18 and C 20 are almost the same, implying that a certain transformation happens at C 18 . The smaller HOP value of C 24 than C 22 is very likely attributed to the finite size effect, because the cationic side chains of C 24 are so long that the current simulation size of 4096 ion pairs is not adequate. Consistent with the RDF results, the HOP of the cationic tail groups increases from C 12 to C 18 , suggesting that the tail groups become more and more aggregated with increasing side-chain length before the transformation. However, the HOP value for the tail groups of C 20 is much smaller than both C 18 and C 22 , while the first COM RDF peak of C 20 is higher than C 18 . This inconsistency may be explained by the fact that the side chains of C 20 are more aggregated after the transformation, while the tail groups are in half way of adjusting from separately aggregated nonpolar domains to layered structures, and thus their aggregation degree is less than both structures represented by C 18 and C 22 , respectively. The decrease of the HOP values for the tail groups of C 24 is again likely due to the finite-size effect. The orientation correlation functions (OCFs) 20 were also calculated for the COMs of cationic side chains to quantify their degree of parallelization. As shown in Fig. 2f, the first peaks of the OCFs for C 12 to C 24 are all located at around 0.5 nm, the distance between two neighboring side chains, so the value of the first OCF peak can quantify the orientation correlation between nearest side chains. The first OCF peak increases with side-chain length from C 12 to C 16 . The value for C 16 is around 0.7, a quite strong orientation correlation, indicating that neighboring side chains are well aligned in parallel. However, there are only two peaks formed in the whole range and the OCF value quickly approaches 0 beyond 2 nm, indicating that no long-range orientation correlations exist for side chains. Therefore, from C 12 to C 16 , with increasing side-chain length, more and more side chains align in parallel locally but no global parallelization is formed. From C 18 to C 24 , the heights of the first OCF peaks are all around 0.9, demonstrating that neighboring side chains are well aligned in parallel. At the same time, multiple peaks form and the OCF values approach a non-zero value with increasing distance, demonstrating that the systems now have a certain degree of long-range order.
To have a better understanding of the degree of global parallelization of side chains, the distribution of the angle between two side chains, whose definition is illustrated in Fig. 3a, have also been calculated and are plotted in Fig. 3b. As shown in Fig. 3a, the direction of a side chain is defined as pointing from its head group (CG site A) to its tail group (CG site E). The angle θ is the twist angle between two side-chain directions. It can be seen that www.nature.com/scientificreports www.nature.com/scientificreports/ from C 12 to C 16 the distribution of θ is quite flat, indicating that the orientation correlation between side chains is very weak. The angle distribution of C 18 has much larger probabilities for smaller angles, and those for C 20 to C 24 are mostly populated in the range of 0-30°, indicating that the side chains are better aligned in parallel for those long chain systems.
Although it is obvious that an abrupt change happens at C 18 for all the above quantities, none of which can unambiguously identify whether there happens a phase transition at C 18 or not. In other words, traditional methods involving correlation functions of spatial symmetries cannot characterize the complex phase behaviors of ILs and ILCs. On the other hand, we will show below that a phase transition at C 18 can be clearly identified by looking at the percolation feature of cationic side chains.

Percolation phase transition.
A cluster is formed by a set of side chains, in which each side chain is "connected" with one or more side chain(s) in the same set. See the Method section for the details of cluster definition. Figure 4a-c show several largest clusters in C 12 , C 16 , and C 22 systems. Figure 4a shows that only a few small clusters form in C 12 since the side chains have weak tendencies of aggregation and parallel alignment. As shown in Fig. 4b, due to growing abilities of aggregation and parallel alignment, more and larger clusters form locally, but the orientations of the clusters are still random. By contrast, for C 22 shown in Fig. 4c, it is apparent that side chains are in parallel globally and the largest cluster is comparable to the system size, while the second largest cluster is very small.
The normalized sizes of the first and second largest clusters with respect to the side-chain length are plotted in Fig. 4d. From C 12 to C 16 , the sizes of these two clusters are both as small as around 1% and slowly grow with increasing side-chain length. The size of the first largest cluster in C 18 grows significantly to be around 25%, and the second largest cluster also increases to be around 5%. For the systems from C 20 to C 24 , the sizes of the largest clusters are all around 70%~80%, comparable to the system size, while the second largest cluster decreases to be around 1%. The slight decrease of the largest cluster size in C 24 is possibly an artifact caused by the limited simulation size. According to the percolation theory 1,3 , the appearance of a giant cluster suggests that a percolation phase transition happens at C 18 . To verify if the abrupt change at C 18 is really a percolation phase transition, we have further calculated both the average cluster size 1 (Fig. 4e) and the correlation length 1 (Fig. 4f) of the side-chain COMs as a function of side-chain length. As shown in Fig. 4e,f, both the average cluster size and the correlation length reach their maxima of around 500 and 0.42, respectively, at C 18 with large fluctuations, indicating that a percolation phase transition does happen at C 18 . Consequently, it is clarified that the NSL to ILC phase transition in ILs is a critical phenomenon, i.e., a second order phase transition.

Discussion
In this work, we have performed the CG MD simulations for the [C n MIM][NO 3 ] (n = 12, 14, …, 24) IL systems with a system size of 4096 ion pairs. To better adapt asymmetric structural feature of ILC, an anisotropic barostat has been applied to all three dimensions, namely, the simulation box size in each dimension is allowed to change independently in the NPT equilibration processes. Our simulation results indicate that, from C 12 to C 16 , separate small clusters composed of cationic side chains aligned in parallel form locally with the average cluster size gradually increases with the side-chain length, while the whole system structure still retains in the NSL state. A percolation happens at C 18 when most of the separate clusters are suddenly unified, corresponding to a large-scale parallel alignment of cationic side chains. The maximization of both the average cluster size and the correlation length as well as their fluctuations manifests that this percolation phenomenon is really a phase transition. The percolation of the cationic side chains in C 20 -C 24 corresponds to the ILC structure with a global long-range correlation. (a) The direction of a side chain is defined as the vector pointing from its CG site A (head) to E (tail), d is the COM distance between two side chains, and θ is the twist angle between two side-chain directions. (b) The ensemble-averaged angle distribution of θ for the systems from C 12 to C 24 . The distribution of θ is quite flat from C 12 to C 16 , suggesting a low degree of parallel alignment. In C 18 , the distribution of θ is much more tilted to smaller angles than the shorter-chain systems. From C 20 to C 22 , the distribution of θ is mostly populated to angles smaller than 30°, indicating that most side chains are well aligned in parallel.
www.nature.com/scientificreports www.nature.com/scientificreports/ Although an abrupt change can be seen at C 18 in RDF, HOP, and OCF plots, only the data analysis method for the percolation phase transition can unambiguously pinpoint the exact phase transition point at C 18 from the NSL state to the ILC state. This critical phenomenon in ILs should have the same universality as the three-dimensional percolation phase transition listed in Table 1. However, it would be tedious to verify this by directly calculating the critical exponents for IL systems on the basis of performing MD simulations with various system sizes.
Our previous CG MD simulation work 20 has found a sharp transformation around n = 14 at T = 400 K from the NSL state to the ILC state for [C n MIM][NO 3 ] when the cationic side chain length increases, which qualitatively agrees with experimental observations 12,16,17,24,25 . However, the exact transition point of n = 14 was questionable since (1) the simulation size of 512 ion pairs is not large enough to eliminate the finite-size effect, and (2) the isotropic barostat may be inappropriate for ILC simulations. The simulations performed in this work with a much larger size of 4096 ion pairs and an anisotropic barostat clearly shows that the structural change around n = 14 becomes mild, and a sharp percolation phase transition appears at n = 18.
Some experiments 26,27 reported a similar phase transition of ILs changing from the NSL state with side chains locally aligned in parallel to the ILC state with side chains globally aligned in parallel. However, the exact transition point may vary according to different conditions, such as temperature, type of ions, etc. The purpose of this work is not providing accurately the transition point, rather it provides a qualitative and general theoretical framework for understanding the phase behavior of ILs changing from the NSL state to the ILC state with increasing cationic side-chain length, which may be very important for the applications of ILs and ILCs in dye-sensitized solar cells 28 , electrofluorescence switches 29 , electrolytes for Li-ion batteries 30 and electrochemical sensors 31 .
With complex molecular structures and interactions, complex fluids frequently exhibit complex phase behaviors beyond the description of simple liquid theories. Therefore, the identification methods for the phase transition of simple liquids, such as the translation and orientation correlation functions, may not be applicable to determining the phase transitions of complex fluids. In this work, we have demonstrated that, as typical complex fluids, the phase transition of IL systems from the NSL state to the ILC state cannot be unambiguously identified by the correlation functions requiring high spatial symmetry due to the complexity of the spatial features of the NSL and ILC phases. By contrast, the data analysis method for percolation easily and clearly identifies the phase transition point, because the percolation phenomenon merely requires a minimal spatial symmetry of connectivity. These results demonstrate the necessity of developing novel liquid theories for complex fluids as well as the importance of percolation theory in describing phase behaviors of complex fluids. www.nature.com/scientificreports www.nature.com/scientificreports/ Methods CG MD simulation. The [C n MIM][NO 3 ] ILs, n = 12, 14, …, 24, are modeled by the EF-CG force field 21,22 . All the simulated systems contain 4096 ion pairs in a parallelepiped box with the PBC applied to all three dimensions. A cutoff distance of 1.4 nm was applied to the VDW and the real part of the electrostatic interactions, and the particle-mesh Ewald method 32 was employed to handle the long-range electrostatic interactions. The temperature was kept constant by using the Nosé-Hoover thermostat 33 with a time constant of 0.5 ps. For the NPT simulated annealing processes from 1200 K to 800 K, the isotropic Parrinello-Rahman barostat 34 with a time constant of 2 ps was employed to keep the system pressure constant. For all other NPT cases, the anisotropic Parrinello-Rahman barostat with P = 1 atm was applied to the three dimensions independently. All the simulations were performed by using the GROMACS software package 35 with a time step of 4 fs. For each IL system, a large random initial configuration was first equilibrated by a 1-ns NPT MD simulation at P = 10 atm and T = 100 K to reduce molecular distances, followed by an NPT simulated annealing process cooling from 1200 K down to 800 K with a temperature interval of 100 K and 2-ns simulated time at each temperature. For convenience, an isotropic barostat was applied for this annealing process. Another subsequent NPT simulated annealing process was performed from 800 K down to 400 K with a temperature interval of 100 K and 8-ns simulated time at each temperature, during which an anisotropic barostat was employed to adjust the size in each dimension independently. The systems finally went through an equilibrium NPT (anisotropic barostat) run at T = 400 K for 8 ns with 2000 configurations evenly sampled from the simulation trajectory.
Heterogeneity order parameter. The HOP 23 continuously quantifying the spatial heterogeneity is defined as where r ij is the distance between site i and site j corrected by the PBC. The normalized distance parameter σ = (V/N) 1/3 with V being the system volume and N being the total number of sites.
Orientation correlation function. The OCF 20 for side chains were calculated to quantify the degree of parallel alignment of side chains. It is defined as the ensemble-averaged orientation correlation between two side chains as a function of their COM distance: where → u r ( ) i is the unit vector of cation i located at → r i pointing from CG site M1 to site E.

Twist angle distribution.
To quantify the degree of global parallelization, the ensemble-averaged twist angle distribution between two side chains were calculated as The twist angle distribution is close to uniform if all side chains are randomly oriented or locally aligned in parallel but globally randomly oriented. Only when side chains are globally aligned in parallel can the angle distribution be mostly populated in small angles.

Definition of cluster.
A cluster is defined as a set of "connected" side chains. The COM distance and the twist angle between two side chains are combined together to identify whether they are connected. The COM RDFs in Fig. 2d indicate that the first shell is within 0.72 nm, and the twist angle distribution in Fig. 3b shows that the angle is mostly within 30° for high parallelization cases. Therefore, two side chains are considered "connected" if the COM distance is less than 0.72 nm and at the same time the twist angle between the two side chains is less than 30°. In a cluster, each side chain is connected to at least one side chain in the same cluster.
Average cluster size and correlation length. The probability that a site belongs to cluster i of size S i is S i /N, so the probability that a randomly chosen cluster has the size of S i is S N / i 2 . Therefore, the average cluster size S ave of a given system can be calculated by 1 where N c is the total number of clusters. Note that theoretically only sites belong to finite clusters should contribute to the calculation since an infinite cluster would make the average cluster size diverge, so numerically the largest cluster after percolation phase transition is not counted in the calculation of the average cluster size. The correlation length ξ is defined as the ensemble average of the distance between two sites in the same finite cluster 1 : www.nature.com/scientificreports www.nature.com/scientificreports/ where the radius of gyration R i of cluster i is defined as = ∑ − R r r S ( ) / i j k j k i 2 1 2 , 2 2 . The correlation length ξ is indeed directly related to the average cluster size S ave , because a longer correlation length always corresponds to larger clusters given that only pairs in the same cluster contribute to the correlation length.