Coarse grained model for calculating the ion mobility of hydrocarbons

Hydrocarbons are widely regarded as insulating compounds in industries. Understanding electrical properties and improving insulation performance may contribute to developing powerful electric applications. However, fundamental characteristics of hydrocarbons in conduction phenomena, such as their ionic behavior, are not completely understood.¹ Under some conditions, convection, turbulence, and chaos are induced. Although hydrocarbons are generally electrically neutral, when ultraviolet light is irradiated on an aluminum cathode surface inserted in hydrocarbons, ionized hydrocarbons can be released.² In such cases, the ion mobility $μI$ is defined with respect to the proportionality factor of the ionic drift velocity $vd$ in electric field $E$ as

The ion mobility is an important parameter in the following equation pertaining to the electric current in dielectric liquids:

where j is the current density, q is the charge density, and D is the self-diffusion coefficient. Thus, we can understand the ionic behavior or current in hydrocarbons by determining the ion mobility. In addition, the ion mobility is known to follow Walden’s rule:³ 

where $η$ is the viscosity and A is a constant value.

Many researchers have measured negative and positive ion mobilities in a variety of liquids.^4–8 They have found that the constant value in Walden’s rule for an singly charged anion is around -3.0×10⁹ A·s/cm and that for a singly charged cation is around 1.5×10⁹ A·s/cm.⁹ These experimental studies provide important fundamental information. However, complicated procedures are required in addition to tight control of the temperature and the purity of the liquids, which drastically influence the behavior of ions.¹⁰ Hence, a great deal of effort is required to determine reasonable ion mobilities from experiments and this is the largest problem to understand the conduction phenomenon in hydrocarbons.

Molecular dynamics (MD) simulation is a method that traces the motion of each molecule. MD simulations can effectively analyze sensitive phenomena affected by their surroundings.^11,12 Furthermore, MD simulations provide detailed information about events at the atomic scale, meaning that the mechanisms of phenomena can be easily understood by giving theoretical interpretations of experimental results through MD simulations.^13,14 In this study, we theoretically calculated the ion mobilities of n-hexane, n-heptane, n-octane, and n-nonane from MD simulations and investigated their characteristics.

The coarse grained (CG) model used in this study assumes a molecule of hydrocarbon to be a bead, as shown in Fig. 1.

The center of mass of the molecule is at the center of the bead. The CG force field is generally defined by nonbonded and bonded interactions.¹⁵ However, we do not calculate bonded interactions because we assume a molecule to be a single bead.

We calculated intermolecular interactions using the Lennard-Jones (LJ) potential:

where $rij$ is the separation, $σij$ is the closest distance between two molecules, and $εij$ is the strength of their interaction. Note that the LJ parameters have not been clarified for all molecules, meaning that we have to determine unknown potential functions and parameters for the target molecules in advance. In our model, we simply defined $εij$ and $σij$ as the sum of the atoms contained in the molecules. We used transferable potentials for phase equilibria (TraPPE)¹⁶ to assess $εij$ and $σij$ of each united atom. The parameters of the force field are shown in Table I.

We used a unit cell of hydrocarbons with a size of $2σ×2σ×2σÅ3$ and a CG representative volume element (RVE) consisting of 256 beads. To reduce the computational time required to reach an equilibrium state, we defined the initial configuration of molecules as a face-centered cubic lattice in which each molecule interacts. However, in the case where the crystal structure is used as the initial configuration placement, the total interaction force is balanced because of its symmetric structure. Hence, initial velocities are required. In an equilibrium state, the velocities follow a normal distribution. Therefore, we defined the initial velocities as random values, which follow a normal distribution by the Box-Muller method:

where $vs$ is the initial velocity, and u₀ and u₁ are random numbers between 0 and 1. Furthermore, to conserve the mass density of the system, we set periodic boundary conditions. The initial configuration of molecules in this study is shown in Fig. 2.

To reduce the effect from thermal variation, we determined the ionic drift velocity $vd$ as the time-averaged ionic drift velocity $⟨vd⟩t$⁠. We then calculated the ion mobility $μI$ as

When we apply an electric force field, the velocity of ions has constant acceleration. However, ionized molecules encounter resistance from surrounding molecules in the opposite direction of the electric field. Thus, the ion velocity converges to a steady value. To take this effect into account, we calculated $μI$ as

where u is the bulk velocity.

In this study, the system was allowed to equilibrate with the NVT ensemble at room temperature of 293.15 K. In the NVT simulations, we used the velocity scaling method¹⁸ and the Nose-Poincare method¹⁹ to control the system temperature. The time step was set to 4 fs. We applied the velocity scaling method for 1 ps and the Nose-Poincare method for 200 ps.

According to Walden’s rule, the product of the ion mobility and the viscosity is constant. However, some previous studies have found that the following equation is more suitable for compounds:¹⁷ 

Equations (3) and (7) show that the viscosity of a liquid is an important parameter to determine the ion mobility. We firstly investigated the dependence of the ion mobility on the viscosity.

We ionized the particle located closest to the center of the system, and determined its ion mobility by applying an electric field to the system. The applied electric field E was 30×10⁴ V/cm. Note that the influence of the ionized particle on the applied electric field E is not calculated, because our system contains only one ionized particle. Furthermore, the size of the ionized particle is assumed to be the same as the neutral particles in this study, because the size variation of a single particle is small compared with the size of the system. The time evolution of the ion mobilities of the hydrocarbons are shown in Figs. 3(a)–3(d). The ion mobilities of the particles start to converge to a certain value at around 80 ps and have converged at 200 ps. Comparison of the experimental data from Walden’s rule, Adamczewski,²⁰ Gzowski,²¹ and the data from our MD simulations is shown in Table II. The calculated ion mobilities fairly well agree with the measured values. The dependence of the ion mobility on the reciprocal of the viscosity determined by each approach is shown in Fig. 4. From this figure, the relation between the ion mobility and the reciprocal of the viscosity can be expressed as

Therefore, the gradient of each line in Fig. 4 is the constant value in Eq. (9). The gradients calculated by the least-square method are shown in Table III. The gradient calculated from the computational results is consistent with those obtained from experiments. It is difficult to identify the cause of the differences of the ion mobilities because both experiments and calculations have uncertainties.

In experiments, contamination of the liquids and the uncertainty of phenomena will cause errors in the measurements. In the simulation, the LJ parameters may have some errors. Furthermore, our results include errors owing to the simple model of the molecules, scale size, and lack of Coulomb interactions. However, our results are within the range of experimental values, which means that these factors may not influence calculation of the ion mobility and understanding the fundamental characteristics of conduction phenomena.

The effect of the electric field on the ion mobility remains the subject of dispute, i.e., some researchers report that the ion mobility must depend on the electric field²² and others report that it must not.²³ Both views are based on experimental data. The fact that opposite conclusions are drawn for the same phenomenon indicates the difficulty of experiments under a high electric field, meaning that MD simulations can play an important role in understanding this phenomenon.

We calculated the ion mobility with electric fields of 10×10⁴, 20×10⁴, and 30×10⁴ V/cm. The ion mobility calculated under each electric field is shown in Fig. 5. The ion mobility slightly varies around a constant value for each hydrocarbon. Thermal fluctuation may cause the slight variation in Fig. 5. Temporal variation of the average temperature for n-hexane is shown in Fig. 6. The average temperature varies around room temperature and we can conclude that the NVT ensemble performed well in the simulation. However, in general, the temperature is not constant. The particle speed refers to the Maxwell-Boltzmann distribution and the homoscedasticity of energy. Thus, variation of the temperature may influence the molecule velocity. The motion of molecules caused by such heat fluctuation is a probabilistic phenomenon. Thus, the influence of this probabilistic phenomenon on the ion mobility is based on the central-limit theorem. The standard error of the ion mobility $SeμI$ resulting from thermal fluctuation can be calculated by

where $sμI$ is the standard deviation of the temperature and N_t is the time step. It is obvious that $sμI$ decreases with increasing N_t. The slight deviation of the calculated ion mobility in Fig. 5 must be caused by this standard error, although it can theoretically be reduced and is negligible when the time step increases. From this discussion, we conclude that the ion mobility does not depend on the electric field.

There has been a lot of discussion about the influence of temperature on chemical reactions.^24,25 In particular, the reaction rate changes with temperature.²⁶ The ion mobility, which causes conduction, may change with temperature because conduction is just a type of chemical reaction such as electrode reaction. Thus, we can estimate errors by the temperature and confirm the appropriateness of the measurement by evaluating the dependency of the ion mobility on the temperature.

We calculated the ion mobilities of the hydrocarbons in the temperature range between 273.15 and 373.15 K with an interval of 20 K. Then, we constructed Arrhenius plots of the calculated ion mobilities at each temperature. We then calculated the activation energy $EμI$ from the Arrhenius equation:

where $μI0$ is the frequency factor and k_B is Boltzmann’s constant. We determined the self-diffusion coefficients D from Einstein’s relation:

where e is the elementary charge. We then calculated the activation energy of diffusion E_D by a similar process to that of $EμI$⁠.

Arrhenius plots of the ion mobilities and self-diffusion coefficients of the four hydrocarbons are shown in Figs. 7(a) and 7(b), respectively. The figures show that the ion mobility and the self-diffusion coefficient exponentially increase with increasing temperature in each liquid. This means that the temperature may cause measurement errors. In other words, the change in the temperature during experiments makes measurements complicated.

In addition, Figs. 7(a) and 7(b) show this reaction mode is Arrhenius type, suggesting that we can estimate the ion mobility at any temperature by

where $μI1$ and $μI2$ are the ion mobilities at temperatures T₁ and T₂, respectively. The ion mobilities and self-diffusion coefficients at 293.15 K and their activation energies in each hydrocarbon are given in Table IV. Table IV shows that the ion mobility and self-diffusion coefficient decrease with increasing molecular radius. This is because of the increase in viscous resistance to solvents with increasing size of the molecule. This tendency agrees with the experimental results of previous studies.^27–29 Additionally, the activation energies of the ion mobility and self-diffusion coefficient increase with increasing molecular size. From the above discussion, we can conclude that the reactivity decreases and the energy of ionization increases with increasing molecular size.

We have theoretically calculated the ion mobilities of hydrocarbons and investigated the effects of viscosity, electric field, and temperature on the ion mobility using the CG model in MD simulations. The results show that the ion mobility follows Walden’s rule in the case of hydrocarbons. In addition, the ion mobility depends on the viscosity, but it does not depend on the applied electric field. The ion mobility and the self-diffusion coefficient exponentially increase with increasing temperature, while their activation energies decrease with increasing molecular size.

The main purpose of this study was to understand the fundamental characteristics of conduction phenomena correctly using the CG model. However, more accurate physical properties and information can be obtained using united atom or all atom models. In the future, we will use these models to better understand not only the general phenomena in conduction, but also to predict various phenomena in electrochemistry such as solution chemistry and surface chemistry.