Surface forces mediated by room-temperature ionic liquids (RTILs) play an essential role in diverse applications including self-assembly, lubrication, and electrochemical energy storage. Therefore, their fundamental understanding is critical. Using molecular simulations, we study the interactions between two nanorods immersed in model RTILs at rod-rod separations where both structural and double layer forces are important. The interaction force between neutral rods oscillates as the two rods approach each other, similar to the classical structural forces. Such oscillatory force originates from the density oscillation of RTILs near each rod and is affected by the packing constraints imposed by the neighboring rods. The oscillation period and decay length of the oscillatory force are mainly dictated by the ion density distribution near isolated nanorods. When charges are introduced on the rods, the interaction force remains short-range and oscillatory, similar to the interactions between planar walls mediated by some protic RTILs reported earlier. Nevertheless, introducing net charges to the rods greatly changes the rod-rod interactions, e.g., by delaying the appearance of the first force trough and increasing the oscillation period and decay length of the interaction force. The oscillation period and decay length of the oscillatory force and free energy are commensurate with those of the space charge density near an isolated, charged rod. The free energy of rod-rod interactions reaches local minima (maxima) at rod-rod separations when the space charges near the two rods interfere constructively (destructively). The insight on the short-range interactions between nanorods in RTILs helps guide the design of novel materials, e.g., ionic composites based on rigid-rod polyanions and RTILs.

## I. INTRODUCTION

The past decades have witnessed extensive interest and attention to room-temperature ionic liquids (RTILs) due to their unique combination of various properties such as low volatility, wide electrochemical windows, and excellent thermal stability.^{1–3} Because of these properties, RTILs have found many potential and even practical applications in diverse areas including electrochemical energy storage, green solvents, and drug delivery, to name a few.^{1,4–6} In many of these applications, RTILs are confined in nanostructures or between extended surfaces, e.g., they are confined between nanorods in macromolecular ionic composites and ion gels^{7} or in nanopores in porous electrodes of supercapacitors^{2,4,8} or between extended surfaces when used as lubricants.^{9–13} Because the structure of the interfacial RTILs and the surface forces mediated by RTILs play an essential role in the aforementioned applications, they have been intensively investigated in recent years.

The interactions between surfaces mediated by RTILs are often investigated experimentally using surface forces apparatus (SFA) and atomic force microscopy (AFM).^{14,15} The forces between two negatively charged mica surfaces separated by ethylammonium nitrate were measured using SFA as early as 1988.^{16} Systematic experimental studies using various surfaces (e.g., mica and silica) and RTILs (e.g., 1-ethyl-3-methylimidazolium ethylsulfate, 1-butyl-3-methylimidazolium bis(trifluoromethanesulfonyl)amide, etc.) have been reported recently.^{11,17–20} These studies revealed that, for small surface separations of typically less than 10 ion diameters, the surface force resembles the classical structural force, i.e., the force is oscillatory with its period commensurate with the estimated diameter of an ion pair and it decays rapidly as the surface separation increases. While some experiments suggest that surface forces mediated by some RTILs, especially protic ones, are short-range, a growing body of experiments shows that such surface force can also be long-range, especially when aprotic RTILs are involved.^{20–23} Specifically, the surface force decays exponentially as the surface separation increases and decay lengths of 5-10 nm have been reported.^{15} Such long-range interactions are not anticipated from the classical DLVO theory. Although their physical origins are not fully understood, it has been suggested that they are electrostatic in nature based on the observation that they depend greatly on the ion density and the dielectric constant of the RTILs.^{15,21,24}

Since the inter-surface forces, especially their short-range oscillations, are closely related to the organization and distribution of RTILs confined between surfaces, the structure of RTILs can be inferred from the force-distance curves measured experimentally. For instance, based on the oscillation of the interaction force, it has been suggested that eight to nine ordered layers of ions form in the RTILs that are confined between mica surfaces,^{16} although a more elaborate assignment of the RTIL structure is also possible.^{25} The layered structure of RTILs near solid surfaces has been confirmed experimentally using high-energy x-ray reflectivity studies.^{26} However, the direct measurement of the RTIL structure, especially in confined space, is usually challenging. This limitation, however, can be addressed partially using theoretical or simulation studies in which the structure of RTILs is resolved explicitly at molecular scales. Indeed, the structure of RTILs confined between surfaces (usually planar walls) has been investigated extensively in the past decade using analytical theories, molecular simulations, and classical DFT calculations.^{27–40} The findings of these studies are largely in reasonable agreement with those derived from experimental studies, and interested readers are referred to two recent reviews for further details.^{2,3}

While the theoretical and computational studies of the structure of confined RTILs are fairly advanced by now, similar studies of the surface forces mediated by RTILs are still at a relatively early stage. The classical DLVO theory is not expected to describe the interactions between charged surfaces separated by RTILs because the excluded volume effects associated with the finite size of ions and the correlations between ions, both important for RTILs, are not taken into account in this theory.^{10} In addition, the Debye length, the characteristics length scale for the electrostatic inter-surface interactions in the DLVO theory, is usually smaller than the ion size in RTILs. This suggests that other length scales, e.g., the ion-ion correlation length,^{31,41} are involved in the long-range electrostatic interactions revealed in recent experiments. The limitations of the classical DLVO theory can be addressed using more advanced theories or molecular simulations, which have been successful in predicting the structure of simple liquids, aqueous electrolytes, and polymers confined between surfaces and inter-surface interactions.^{42–45} Nevertheless, only a limited number of analytical theories and simulations have been reported for inter-surface interactions mediated by RTILs. Using the Coulomb gas model, Lee *et al.* showed that inter-surface force between two planar walls separated by RTILs decays in an oscillatory manner as the inter-surface separation increases and the force maxima correspond to the insertion of new ion layers into the gap formed between the planar walls, which agree qualitatively with available experimental data.^{46} A key insight from the model is that the inter-surface force is controlled strongly by the fugacity of bulk RTILs. More recently, using models built upon the Ginzburg-Landau theory and molecular simulations, Vaikuntanathan, Talapin, and their colleagues studied the interactions of planar surfaces and nanocrystals separated by a high-temperature molten inorganic salt KCl.^{47} Their analytical model and simulations captured the oscillatory nature of the inter-surface forces well. A key insight from their study is that the interference between the charge density oscillations near opposing surfaces greatly affects the sign and magnitude of the inter-surface interactions.^{47} The fact that theories and simulations can capture key aspects of the inter-surface forces mediated by ionic liquids revealed in the experiments suggests that theories and simulations are promising tools for studying these forces.

The prior studies on the interfacial structure of RTILs and surface forces mediated by RTILs have greatly advanced our understanding of these problems. Nevertheless, some important knowledge gaps still remain. For instance, few simulations have been performed to delineate the surface forces and interfacial RTIL structure concurrently to provide insight into their correlation. Furthermore, most of the available studies deal with extended surfaces like planar walls. Situations in which the critical dimension of the objects is comparable to the ion size, which are often encountered in molecular self-assembly, are rarely explored. It is not yet clear to what extent the interactions between these objects resemble the interactions between extended surfaces and therefore many questions remain open. For example, how does surface charge affect the interactions between these objects? What is the range of these interactions? How do the distributions of ions around the objects determine these interactions? In this work, we investigate the interactions between nanorods mediated by RTILs using MD simulations to address these questions. We focus on the situation in which nanorods are fully aligned because of its relevance to the recently synthesized macromolecular ionic composites, in which polyanion rods are highly aligned.^{7} The insight gained here will help lay a foundation for understanding the interactions between tilted rods and surfaces.^{47,48} The rest of the manuscript is organized as follows. The simulation system, molecular model, and methods are presented in Sec. II. The simulation results on the interactions between neutral or charged nanorods and their relation with the underlying RTILs structure near the nanorods are presented in Sec. III. Finally, conclusions are drawn in Sec. IV.

## II. SIMULATION SYSTEM, MOLECULAR MODELS, AND METHODS

Figure 1 shows a schematic of the simulation system, which features two aligned nanorods immersed in model RTILs. The rod-rod distance *D* is defined as the distance measured between the axes of the two rods. We performed two sets of studies, one with neutral nanorods and the other with charged nanorods. For each set of study, the distance *D* was varied systematically to study the rod-rod interactions mediated by the RTILs. The number of ions inside the system was fixed in all simulations. Due to the large size of the simulation box (see Fig. 1 for the dimension of the simulation system), the nanorods occupy less than 0.9% of its total volume. Consequently, a bulk like behavior is recovered at positions away from the rods and ions in the system have approximately the same chemical potential in systems with different rod-rod spacings.

In the spirit of the coarse-grained models used in many prior studies of RTILs, which have been successful in capturing key features of the structure and dynamics of bulk and interfacial ions,^{32,33,36,49} herein we adopted similar coarse-grained models for both the nanorod and model RTILs. Each rigid nanorod consists of 33 layers of carbon atoms distributed evenly in the rod length direction. Each layer includes seven carbon atoms arranged as a regular heptagon with a side length of 0.17 nm. For charged nanorods, small partial charges are assigned equally to each atom in the nanorod to give a line charge density of −2.38 *e*/nm. The structure and charge of the rod thus constructed are relevant to the rigid-rod polyanion poly(2,2′-disulfonyl-4,4′-benzidine terephthalamide), or PBDT, used in the synthesis of some macromolecular ionic composites.^{7} Cations and anions in our model RTILs, both rigid, have geometries identical to the BF_{4}^{−} ions (see Fig. 1) and only their central atoms carry a unit positive/negative charge. To account for ions’ electronic polarizability, a background dielectric constant of *ϵ*_{r} = 2.5 was used in the calculation of electrostatic interactions. All atoms in the nanorod have Lennard-Jones (LJ) parameters of *σ*_{LJ} = 0.355 nm and *ϵ*_{LJ} = 0.293 kJ/mol. Atoms in the cation and anion have LJ parameters of *σ*_{LJ} = 0.336 nm and *ϵ*_{LJ} = 0.360 kJ/mol. More details of the force field for RTILs can be found in Table S1 in the supplementary material.

With the above force fields, the bulk density of the RTILs is *ρ*_{+} = *ρ*_{−} = 4.11 nm^{−3} (P = 1 atm; T = 400 K). The radius of the cation/anion is ∼0.22 nm based on the radial distribution function of bulk RTILs (see Fig. S1 in the supplementary material). The effective radius of each nanorod is ∼0.35 nm based on the ion density profile around isolated nanorods (see Fig. S2 in the supplementary material). The coarse-grained models adopted here do not take into account many of the subtle details of real polyanions and RTILs, e.g., the discreteness of polyanions’ surface charge groups and the complex shape of the ions. With these simplifications, we seek to reveal the most essential features of the interactions between nanorods mediated by RTILs, unobscured by the chemical complexity of RTILs and nanorods. The potential limitations of the coarse-grained model are discussed in Sec. IV.

Simulations were performed using the Gromacs code.^{50} All simulations were executed in the NVT ensemble (T = 400 K). Similar to many prior studies on RTILs, an elevated temperature was used to ensure that the phase space is explored effectively within the time scale accessible to our simulations (∼50 ns). The temperature of the system was maintained using the velocity rescaling thermostat.^{51} The non-electrostatic interactions were computed by direct summation with a cut-off length of 1.5 nm. The electrostatic interactions were computed using the Particle Mesh Ewald (PME) method. The real space cutoff and FFT spacing were 1.5 and 0.13 nm, respectively. For each system, an equilibrium run of 20 ns was performed, followed by a 30 ns production run.

During the production run, the force acting on each atom of the rod was recorded and summed to give the total force acting on each rod. The free energy (potential of mean force, or PMF) for rod-rod interactions was calculated using $ED=\u2212\u222bD0Df(x)dx$, where *f*(*x*) is the force at a rod-rod separation of *x*, *D*_{0} = 4 nm is the largest rod-rod separation explored in our simulations, and the free energy at rod-rod separation of *D*_{0} is taken as zero. Hereafter, energy is measured in *k*_{B}*T* with *T* = 400 K for convenience, where *k*_{B} is the Boltzmann constant.

## III. RESULTS AND DISCUSSION

### A. Interactions between neutral nanorods

Figure 2(a) shows the rod-rod interaction force as a function of the distance between two neutral nanorods. As the two rods initially in physical contact are separated from each other, they first attract each other (0.7 nm < *D* < 1.0 nm). As the separation further increases, their interaction force exhibits strongly damped oscillations, which nearly vanish at separations larger than ∼2.0 nm. The force profile *f*(*D*) at *D* > 0.9 nm can be fitted to a damped sine wave function

where $Ldf$ is the decay length, $Lpf$ is the oscillation period, and *a* and *b* are fitting constants. Using the data shown in Fig. 2(a), $Ldf$ and $Lpf$ are determined to be 0.50 nm and 0.48 nm, respectively.

The integration of the force shown in Fig. 2(a) leads to the interaction free energy *E*(*D*) between the nanorods presented in Fig. 2(b). *E*(*D*) is also oscillatory and decays with *D* rapidly. Its decay length ($LdE=0.44\u2009nm$) and oscillation period ($LpE=0.51$ nm) are determined by fitting *E*(*D*) to the form of Eq. (1), and they are found to be comparable to those of the interaction force. It can be seen that separating two nanorods initially in contact must overcome an energy barrier of ∼10 *k*_{B}*T*/nm. Such a large energy barrier is attributed mostly to the attractive rod-rod interactions at separations of 0.7-1.0 nm, which are caused by the van der Waals attractions between the two rods and the capillary effects. For *D* ≤ 1.0 nm, one cannot even fit one ion between the closest surface points of the two rods (see Fig. S3 in the supplementary material). As such, a vacuum is developed between the two rods and a vacuum-RTIL interface is formed. To lower the free energy associated with the vacuum-RTIL interface, an attractive force between the rods, which pulls the rods closer to reduce the area of the vacuum-RTIL interface, is developed.

The short-range, oscillatory forces between two neutral rods observed here resemble the structural forces between extended surfaces separated by simple fluids. The latter is usually understood using the contact value theorem and the superposition approximation method.^{10,14} Following the contact value theorem, the interaction force *P* between two planar walls separated by a distance *D* is given by *P*(*D*) = *k*_{B}*T*[*ρ*_{D}(0) − *ρ*_{∞}(0)], where *ρ*_{D}(0) is the contact fluid density (i.e., fluid density on a wall’s surface) when the walls are separated by a distance *D*. Therefore, repulsive (attractive) forces are accompanied by the highly (weakly) ordered fluid structure between the walls and accordingly the high (low) fluid contact density on the walls.^{14} The fluid structure between the two walls, embodied in the contact fluid density and affected by the geometrical constraints imposed by the two walls, can be qualitatively inferred from the fluid structure near isolated walls by superposition:^{52} the fluid density between two walls separated by a distance of *h* varies as *ρ*_{h}(*z*) ∼ *ρ*_{∞}(*z*) + *ρ*_{∞}(*h* − *z*), where *z* is the distance from one wall and *ρ*_{∞} is the fluid density near an isolated wall. It follows that $PD\u223ckBT\rho \u221eD\u2212\rho b\u223ckBT\rho \u221eex(D)$, where *ρ*_{b} is the density of bulk fluids and $\rho \u221eexD=\rho \u221eD\u2212\rho b$ is the excess fluid density at a distance of *D* from the wall. This relation has been used successfully to understand structural forces.^{43,53} Importantly, it suggests that the oscillation period of structural forces between two planar walls is close to the oscillation period of the fluid density near an isolated planar wall, which in turn is comparable to the size of the fluid molecules. The structural forces between curved surfaces are often understood from those between planar walls with the help of the Derjaguin approximation, and it can be shown that key characteristics of structural forces, e.g., their oscillation period, are not affected by the curvature of the surfaces.^{14}

For the interactions between molecularly thin, neutral nanorods in RTILs considered here, the Derjaguin approximation is no longer rigorous since the radius of nanorods, the size of ions, and rod-rod separations are all comparable. Nevertheless, we shall examine to what extent the rod-rod interactions can be understood using the above framework established for the structural forces between extended surfaces. Assuming that the rod-rod interaction is controlled primarily by the structure of RTILs along the centerline connecting the two rods’ axes, we have *f*(*D*) ∼ (*ρ*_{∞}(*D*) − *ρ*_{b}). Because the cations and anions are symmetric in our RTILs and the nanorods are neutral, the fluid densities *ρ*_{∞}(*D*) and *ρ*_{b} are taken as the summation of cation and anion density. Figure 3 compares the rod-rod interaction force *f*(*D*) and the excess fluid density $\rho \u221eexD$. Excellent commensuration between *f*(*D*) and $\rho \u221eexD$ is observed. Importantly, the oscillation period of the rod-rod interactions (0.48 nm) agrees with that of the fluid density near isolated rods (0.48 nm), and both oscillation periods are close to the ion diameter (∼0.44 nm). These results indicate that the structural forces between molecularly thin nanorods immersed in RTILs can be understood qualitatively using the established method for the structural forces between planar walls mediated by simple fluids.

The above discussion suggests that the liquid structure around nanorods, which is controlled by the geometrical constraints imposed by the nanorods, greatly affects the rod-rod interactions. Nevertheless, the superposition method for inferring the liquid structure near the rods does not take into account the rods’ small radius and thus cannot provide a full picture of the liquid structure near them. To address this limitation, next we examine the distribution of RTILs around the nanorods at rod-rod separations corresponding to the three extrema of *f*(*D*) in Fig. 2(a).

The left column of Fig. 4 shows the fluid density distribution (*ρ*_{f} = *ρ*_{+} + *ρ*_{−}) near the two nanorods corresponding to points A, B, and C in Fig. 2(a). Clearly, each nanorod affects the fluid structure near the other rod, especially in the gap between them. The schematics in the right column of Fig. 4 highlight the evolution of the fluid structure as the rods are gradually separated from each other and the areas in which strong interference of fluid structure near neighboring nanorods is enclosed by the red dashed lines. At a rod-rod separation of 1.10 nm, the strong interference of the first ion density peaks centering on the two rods leads to greatly enhanced ion density along a portion of the rod surface [see the high density spots in Fig. 4(a1) and the magenta lines in Fig. 4(a2)], and consequently a strong repulsion between the rods as shown in Fig. 2(a). Figure 4(b2) shows that as the rod-rod separation is increased to 1.30 nm, two kinds of interferences set in: a constructive interference between the first ion density peak centering on each rod and a destructive interference between the first ion density peak centering on one rod and the first ion density valley centering on the other rod. The constructive interference enhances the ion density in the space between the two rods [see Fig. 4(b1) and the grey patch between the two rods in Fig. 4(b2)] but the destructive interference decreases the ion density along a portion of each rod’s surface [see the magenta lines in Fig. 4(b2)]. Following the contact value theorem, the latter interference reduces the repulsion between the rods, which is consistent with the force valley at point B in Fig. 2(a). As the rods move further apart to a rod-rod separation of 1.6 nm, the geometrical constraints provided by one rod to the ions near the other rod become quite weak. However, a weak constructive interference between the first ion density peak centering on one rod and the second ion density peak centering on the other rod leads to a slight enhancement of the ion density along a portion of each rod’s surface [see the magenta lines in Fig. 4(c2)], and consequently a peak in the rod-rod repulsion force in Fig. 2(a).

### B. Interactions between charged nanorods

Next we turn to the rod-rod interactions when both rods are negatively charged. Here, the double layer forces are expected to become important. Figure 5(a) shows the rod-rod interaction force as a function of the distance between two charged nanorods. The interaction force between charged rods shows damped oscillations similar to those between neutral rods. However, some differences are also notable. First, compared to that between neutral rods, the interactions between charged rods have a longer range (∼3 nm) and the period of force oscillation is larger. Indeed, fitting the force curve in Fig. 5(a) to Eq. (1) gives a period of $Lpf=0.58\u2009nm$ and a decay length of $Ldf=0.80$ nm, compared to $Lpf=0.48$ nm and $Ldf=0.50$ nm for the neutral rods. Second, the first valley of the force curve appears at a rod-rod separation ∼0.30 nm farther compared to that for neutral rods and it is shallower than the first valley in the force curve for neutral rods. This observation is caused by the electrostatic repulsion between the two negative charged rods. Finally, a kink is observed in the force curve at a rod-rod separation of 1.15 nm (point A in Fig. 5). This observation can be attributed to the structural force. At a rod-rod separation of ∼1.15 nm, the first counter-ion layer around both rods merges in the gap between them to enhance the ion density there (see Fig. S4 in the supplementary material). The high density of this counter-ion layer helps reduce the electrostatic repulsion between the rods. Meanwhile, since these ions also serve as the liquid medium for the two interacting rods, their high local density in the gap between the rods contributes to a repulsive force peak. This force peak is similar to the repulsive, structural (entropic) force peak for the interactions between two neutral rods (cf. the force peak at point A in Fig. 2), which occurs at essentially the same rod-rod separation.

Figure 5(b) shows the interaction free energy between the two charged rods as a function of the rod-rod separation. The free energy curve shows oscillation period ($LpE=0.59$ nm) and decay length ($LdE=0.87$ nm) similar to those for the force curve in Fig. 5(a). Despite that both rods carry the same charge, a primary free energy minimum occurs at a rod-rod separation of 1.00 nm. Separating two rods initially positioned at this separation must overcome an energy barrier of ∼4.6 *k*_{B}*T*/nm. Such a barrier can help stabilize nanorods in RTILs and may contribute to the mechanical strength of the polyanion nanorod-based ionic composites developed recently.^{7}

Since the current understanding of the interactions between charged rods mediated by RTILs is quite limited, hereinafter we focus on understanding how the interactions oscillate and decay with the rod-rod separation. Because the decay length and oscillation period are similar for the interaction force and free energy and it is customary to focus on the free energy in the analysis of colloidal interactions,^{14} below we examined the interaction free energy between charged rods. In particular, we focused on how the structure of RTILs near two rods evolves as their separation changes and how the structure change affects the oscillation and decay of the free energy with rod separation. Given that the interactions between charged colloids are most closely tied to the distribution of the space charge around the objects, the structure of the RTILs near the rods will be characterized using the local space charge *ρ*_{e} = *e*(*ρ*_{+} − *ρ*_{−}), where *ρ*_{±} is the local cation/anion number density.

The left column of Fig. 6 shows the space charge distribution near the two rods at three rod-rod separations corresponding to the three local free energy maxima in Fig. 5(b) (points B, D, and F). The interference of the space charge distribution near the two rods is evident in the gap between the two rods. Prior studies have shown that, for space charges enclosed between two charged planes, their distribution can be understood from the simple superposition of the space charge near an isolated plane.^{47} Here we check whether this simple superposition is also reasonable for the space charge near rods whose radius and separation from each other are both comparable to the ion size. To this end, we computed the space charge along the centerline connecting the axes of the two rods by superposing the space charge near isolated nanorods and the results are compared to that measured directly from MD simulations in the right column of Fig. 6. For all three rod-rod separations, the space charge along the rod-rod centerline is captured quite well by the superposition method. A key observation of Figs. 6(a2), 6(b2), and 6(c2) is that the interference of the space charge from the opposing rods is “destructive” in nature, i.e., a space charge peak (valley) associated with one rod is cancelled by a space charge valley (peak) associated with the other rod. Such a destructive interference of the space charges near two rods reduces the positive space charge along the centerline of the two negatively charged rods, thus causing less effective screening of the rods’ negative surface charges by the counter-charges in the gap between them. The weakened screening of the electrostatic repulsion between the like-charged rods drives up the overall interaction energy between them energetically, which helps us to explain why free energy is maximized locally at points B, D, and F in Fig. 5(b).

The interference of the space charge near two rods with separations corresponding to the local minima of the free energy [points C, E, and G in Fig. 5(b)] is shown in Fig. 7. A key difference from the interference shown in Fig. 6 is that the interference is “constructive” in nature here, i.e., a space charge peak (valley) associated with one rod is enhanced by another space charge peak (valley) near the other rod. Such a constructive interference tends to increase the positive space charge density in the gap between the two rods and thus leads to more effective screening of the surface charge of the rods. Consequently, the free energy is minimized locally at points C, E, and G in Fig. 5(b).

In the above discussion (in particular the right columns of Figs. 6 and 7), only the interference of the space charge along the centerline of the two rods is explicitly addressed. Consequently, the effect of the finite radius of the rods on space charge interference is not considered. To examine the interference more accurately, we computed the net space charge in the gap between the two rods for rod-rod separations corresponding to points B to G in Fig. 5. The gap is defined as a rectangle box with a height of 1.0 nm and its two sides are bounded by the axes of the two rods (see Figs. 6 and 7). Table I shows that the net space charge in the gap at points B, D, and F is much smaller than at points C, E, and G. This is in good agreement with the difference expected from the destructive interference and constructive interference for these two sets of points.

Maximal free energy positions . | Minimal free energy positions . | ||||
---|---|---|---|---|---|

Point B | Point D | Point F | Point C | Point E | Point G |

+1.24 e/nm | +1.37 e/nm | +1.71 e/nm | +2.61 e/nm | +2.29 e/nm | +2.12 e/nm |

Maximal free energy positions . | Minimal free energy positions . | ||||
---|---|---|---|---|---|

Point B | Point D | Point F | Point C | Point E | Point G |

+1.24 e/nm | +1.37 e/nm | +1.71 e/nm | +2.61 e/nm | +2.29 e/nm | +2.12 e/nm |

The above results suggest that the interference of space charge near molecularly thin rods can be understood reasonably well by superposition of the space charge near isolated rods. Together with the understanding that the destructive (constructive) interference of the space charge around two neighboring rods tends to increase (decrease) their interaction free energy, it is reasonable to conclude that the interactions between thin rods in RTILs can be qualitatively understood from the superposition of the space charge around isolated rods. In a recent study of the interactions between planar surfaces mediated by molten salts, it has been analytically shown that using the superposition method for evaluating the space charge between surfaces, the oscillation period and decay length of the interactions between two surfaces are the same as those for the space charge profile near the isolated surfaces.^{47} Here we checked whether the same would be true for the thin rods considered in this work. Figure 8 shows the space charge density profile near an isolated rod, where the damped oscillation of the space charge is evident. Fitting the space charge curve to the form given by Eq. (1) gives an oscillation period of $Lpe=0.61$ nm and a decay length of $Lde=0.80$ nm. The computed $Lpe$ agrees with the free energy oscillation period of $LpE$ = 0.59 nm extracted from Fig. 5(b). Physically, the oscillation period of the space charge profile near an isolated rod is similar to the thickness of a counter-ion and co-ion layer positioned next to each other shown in Fig. S5 in the supplementary material. Since earlier experiments showed that each period of interaction force (and consequently interaction energy) oscillation corresponds to a neighboring pair of counter- and co-ion layer being removed between two interacting surfaces,^{10} the agreement of $Lpe$ with $LpE$ is physically sound. The computed decay length of the space charge near isolated rods $Lde$ agrees with the rod-rod free energy decay length of $LdE$ = 0.87 nm extracted from Fig. 5(b). These results thus support the method of inferring rod-rod interactions from the space charge profiles near isolated rods.

The results in Figs. 2 and 5 show that the wavelength (*L*_{p}) and decay length (*L*_{d}) of the damped, oscillatory rod-rod interactions mediated by RTILs depend apparently on the surface charge of the rods. Fundamentally, however, *L*_{p} and *L*_{d} are mainly the properties of the RTILs,^{15,54} and to a lesser extent, they are also affected by the geometry of the interacting surfaces. Two important types of length scales exist in RTILs. One is the size of the ions and the other is *L*_{p} and *L*_{d} of the space charge density wave within them (for bulk RTILs, a charge density wave is observed near individual ions; for RTILs near electrified surfaces, a charge density wave is also observed as one moves from the surfaces toward bulk liquids). *L*_{p} and *L*_{d} depend on the size and the interactions between ions and are thus primarily the properties of the RTILs. For RTILs near electrified surfaces, *L*_{p} and *L*_{d} are also affected by the geometry of the interacting surfaces to some extent because the distribution of space charge near an electrified surface is known to be affected somewhat by its curvature.^{55} For example, *L*_{d} is 1.25 nm in our bulk RTILs, but it is reduced by ∼20% in the same RTILs near planar walls (surface charge density: 0.75 *e*/nm^{2}) and is reduced further by ∼15% in the same RTILs near the charged rods considered in this work.

The principal role of the surface charge is to control which length scales in RTILs are manifested more prominently in rod-rod interactions. The interactions between neutral rods are dominated by structural forces, which are mainly controlled by the packing of ions in the gaps between the rods. Consequently, the wavelength and decay length of these interactions are controlled mostly by the ion size. In comparison, the interactions between highly charged rods are dominated by double layer forces, which depend strongly on the interferences of the space charge near individual rods. As such, the wavelength and decay length of these interactions are controlled primarily by *L*_{p} and *L*_{d} of the space charge density waves within RTILs. One can expect the wavelength and decay length of the rod-rod interaction force/energy to evolve continuously as the surface charge on the rods deviates from zero. A systematic investigation of such an evolution is beyond the scope of the present work but is worth pursuing in the future.

The rod-rod interactions revealed in our simulations are rather short-range. This seems to be in agreement with the range of protic RTIL-mediated surface forces measured experimentally,^{12,15} but different from the long-range interactions reported in many other experiments (e.g., a decay length of ∼5-10 nm has been reported).^{15,54} The absence of long-range rod-rod interactions in our simulations can potentially be caused by the limitations of the MD simulations and may also be related to the nature of the model RTILs and the nanorods used here. First, in our simulations, the periodic boundary conditions are used. Therefore, long-range rod-rod interactions may be smeared out because of the interactions between the rods and their images. Second, since only very strong forces can be reliably measured in MD simulations, it is possible that the weak long-range interactions are simply buried in statistical noises. Finally, some recent theoretical and experimental studies suggested that the decay length of the inter-surface interactions in RTILs is likely a property of bulk RTILs.^{54} Given that the model RTILs adopted here have a small ion diameter (∼0.44 nm) and the nanorods are molecularly narrow, it is possible that inter-surface interactions mediated by our model RTILs are truly short-range. In light of the uncertainties associated with these considerations, we caution that the present study of rod-rod interactions mediated by RTILs does not provide a definitive answer on the decay length scale of RTIL-mediated interactions. Instead, it mainly provides insight into the interactions at short-range when both structural and double layer forces are strong, which are important in applications such as synthesis of macromolecular ionic composites using polyanions and RTILs.^{7}

## IV. CONCLUSIONS

Using molecular simulations, we have studied the interactions between nanorods mediated by RTILs when both structural and double layer forces are important. When the rods are neutral, the force between them shows damped oscillation as a function of increasing separation, and both the oscillation period and decay length of the interaction force are close to the ion diameter. Such oscillatory force is essentially the classical structural force and is controlled by the packing constraints imposed by the rods to their neighboring RTILs. When the rods are charged, the force between them still shows damped oscillation as their separation increases. However, both the oscillation period and the decay length are larger than those for the interactions between neutral rods. In addition, the first valley of their interaction force becomes shallower and shifts to a larger rod-rod separation. Nevertheless, separating two rods initially positioned at the first force valley must overcome an energy barrier of ∼5 *k*_{B}*T*/nm.

We have shown that, using the simple superposition method for estimating the ion density profile between rods and the classical contact value theorem, the oscillation pattern and decay length for interactions between neutral rods can be estimated using the ion density profile despite the diameter of the rod being comparable to the ion size. Likewise, the interactions between charged rods can be interpreted qualitatively by considering the interference of the space charge around different rods. It is found that, even for molecularly thin rods, the interference of space charge can be estimated from the space charge near isolated rods using the simple superposition method. The destructive (constructive) interference of the space charge near neighboring rods leads to the reduced (enhanced) screening of the surface charge of the rods and thus higher (lower) rod-rod interaction energy. Because of such interferences, the oscillation period of the rod-rod interaction energy is close to the thickness of a pair of counter- and co-ion layers, in good agreement with that found experimentally. The rod-rod interactions revealed by our simulations are short-range. However, for reasons such as the nature of the model RTILs used and noise in the measurement of interaction forces, we caution that this result does not preclude the existence of long-range force between charged surfaces separated by RTILs.

We note that the insight gained in this work is derived from simulations based on RTILs featuring cations and anions with identical size. In practice, many RTILs have cations and anions with different sizes (e.g., [BMIM][Cl]). Nevertheless, the basic features of the rod-rod interactions mediated by these asymmetric RTILs can be inferred from the present study. Specifically, the periodicity and decay length (*L*_{p} and *L*_{d}) of the interactions between neutral rods should be determined primarily by the size of the larger ions. This is because the distribution of ions near the rods, which governs the structural force between them, is constrained mostly by the packing of the larger ions. On the other hand, *L*_{p} and *L*_{d} of the interactions between highly charged rods should still depend mostly on *L*_{p} and *L*_{d} of the charge density wave in RTILs. Since *L*_{p} and *L*_{d} of the charge density wave in RTILs are affected greatly by both the small and the large ions, we expect the small ions to have a stronger effect on *L*_{p} and *L*_{d} of the interactions between highly charged rods than between neutral rods. A systematic study of how asymmetric ions affect RTIL-mediated interactions between nanorods will be pursued in the future.

## SUPPLEMENTARY MATERIAL

See supplementary material for the force fields of nanorod and RTILs, the radial distribution function of ions in bulk RTILs, the ion distribution near isolated neutral rod, the ion distribution near neutral rods with a separation of 1.0 nm, the ion distribution near charged rods corresponding to point A in Fig. 5, and ion distribution near isolated charged rod.

## ACKNOWLEDGMENTS

We thank the ARC at Virginia Tech for generous allocations of computer time on the BlueRidge and NewRiver cluster. R.Q. gratefully acknowledges the support from NSF (CBET-1461842). R.Q. was partially supported by an appointment to the HERE program for faculty at the Oak Ridge National Laboratory (ORNL) administered by Oak Ridge Institute of Science and Education. J.H. and B.G.S. acknowledge work done at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. We thank Dr. Alpha Lee for insightful comments.