Molecular Dynamics Simulation of the Interfacial Structure and Differential Capacitance of [BMI+][PF6−] Ionic Liquids on MoS2 Electrode

Abstract

MoS₂ nanomaterials and ionic liquids (ILs) have attracted tremendous interest as the prime electrodes and electrolytes of supercapacitors. However, the corresponding charge storage mechanisms have not yet been clearly understood. Herein, we study the molecular-level energy storage mechanisms of the MoS₂ electrode in imidazolium ionic liquid ([BMI⁺][PF₆⁻]) using molecular dynamics (MD) simulation. The electric double-layer (EDL) structures of MoS₂ electrodes in [BMI⁺][PF₆⁻] electrolytes are comprehensively studied in terms of number density, MD snapshots, separation coefficient, and ion-electrode interaction energy. Based on this, the electric potential and electric field distributions are calculated by integrating Poisson equations. Importantly, a bell-shaped differential capacitance profile is proposed, different from the U-shaped curve from the conventional Gouy–Chapman theory. Especially, it can be well reproduced by the differential charge density curve in the Helmholtz layer. This indicates that the capacitive behaviors of the MoS₂ electrode are primarily determined by the counterion population/structure in the EDL region 5.0 Å from the electrode surface. In the end, ion diffusion coefficients within different interfacial EDL regions are evaluated, revealing that dynamics are significantly suppressed by ~50% in the Helmholtz region. However, the dynamics can be recovered to a bulk state with the ion position 10 Å away from the electrode surface. The as-obtained insights reveal the charge storage mechanisms of MoS₂ in ILs, which can provide useful guidance on improving the energy density of MoS₂ supercapacitors.

1. Introduction

Supercapacitors have attracted great attention due to their high charging/discharging rate, superior power density, and ultralong cycling life, which are widely used in power-needed areas (e.g., uninterruptable demands in wireless sensors, electronic modules, and subway and transit networks) [1,2,3]. However, the limited energy density significantly impedes their applications. The energy density is proportional to the specific capacitance of electrode materials and the square of the electrochemical stability window of the electrolyte. Therefore, developing advanced electrodes composed of nanomaterials with high surface area, excellent electrical conductivity, etc., as well as electrolytes with high electric windows, are the key strategy to improve the energy density of supercapacitors [4,5,6].

Two-dimensional (2D) molybdenum disulfide (MoS₂) has attracted broad interest in energy storage and conversion applications. Among various kinds of MoS₂ structures, 1T phase MoS₂ has been widely investigated as the electrode material of supercapacitors [7,8]. 1T phase MoS₂ is composed of a layer of Mo atoms sandwiched between two layers of S atoms, which is used as the simulation electrode model in this work. Compared to the 2H phase, 1T phase MoS₂ has two unique features that facilitate both electron transport and ion diffusion for electrochemical energy storage technologies. 1T phase MoS₂ has distorted octahedral coordination, enabling the electrical conductivity to be 10⁵ times higher than the 2H phase; meanwhile, its hydrophilic property benefits ion diffusion [9]. For example, M. Chhowalla et al. suggested that 1T phase MoS₂ can achieve capacitance values ranging from ∼400 to ∼700 F cm⁻³ in a variety of aqueous electrolytes [10], which is much higher than that of carbon-based supercapacitors (e.g., graphene and activated carbon, typically 100~300 F cm⁻³). Moreover, they further demonstrated that MoS₂ is also suitable for high-voltage operation in non-aqueous organic electrolytes, showing prime energy and power density values, superior coulombic efficiencies (>95%), and stability (>5000 cycles). Feng et al. reported that 3D conductive MoS₂/N-doped carbon heteroaerogels exhibit ultrahigh capacitance (~4000 F/g), remarkable rate capability (~2000 F/g at 10 A/g), and excellent cycling stability (~8% capacitance decay after 8000 cycles) [11]. Recently, our work realized a high capacitance of 511 F cm⁻³ at an industrial-level submillimeter MoS₂ film electrode (94.2 μm) [12]. These experimental studies have demonstrated the great potential of MoS₂ as the electrode material of supercapacitors.

On the other hand, room-temperature ionic liquids (ILs), composed of bulky organic cations and versatile anions, have been considered as the promising electrolyte of supercapacitors due to the large electrochemical stability window (3~7 V), negligible vapor pressure, superior thermal and chemical stability compared to aqueous and organic electrolytes [13,14]. The capacitive behaviors of ILs are primarily influenced by their electric double-layer (EDL) structure and dynamic properties, which have been widely investigated on graphene electrodes (including ours) [15,16,17,18]. For example, using molecular dynamics (MD) simulations, J. G. McDaniel studied the differential capacitance profiles of ionic liquids ([BMIM⁺][BF₄⁻] and [BMIM⁺][TFSI⁻]) mixed with acetonitrile and 1,2-dichloroethane on graphene electrodes. It is reported that both pure and 10% mole fraction ionic liquid electrolytes exhibit camel-shaped capacitance profiles with two peaks on either side of a minimum at the potential of zero charge (PZC) [19]. Using MD simulation, J. Seebeck reported that the differential capacitance of [BMI⁺][PF₆⁻] shows a fluctuating behavior with respect to the curvature radius of the electrode [18]. E. E. Fileti studied the effect of the ionic size of ILs on the double-layer capacitance [20]. G. Feng investigated the effects of hydrophobicity/hydrophilicity of ILs and carbon electrodes on the interfacial distribution of ions and electrosorbed water [21]. It is revealed that using hydrophilic ionic liquids would help to keep water molecules away from the negatively charged electrodes, even at large electrode polarizations. Utilizing constant voltage MD simulation, Tu et al. suggested that both pure and 10% mole fraction ionic liquid electrolytes exhibit camel-shaped capacitance profiles with two peaks on either side of a minimum centered at the potential of zero charge [19]. More interestingly, Tu et al. investigated the relationship between capacitance and interfacial structure, revealing that inner layer capacitances for the organic solvents and [BMIm⁺][BF₄⁻] solutions are very similar, regardless of solvent type and largely independent of ion concentration [22].

Based on the above analysis, combinations of MoS₂ electrodes and high-stability ILs show great potential to boost the energy density of supercapacitors. However, the underlying charge storage mechanisms are poorly understood at the molecular scale, which needs a comprehensive MD simulation study. On the graphene electrode, conventional U-shaped (from Gouy–Chapman theory), camel-shaped, and bell-shaped differential capacitance curves have been predicted by MD simulation [23,24,25]. It is interesting to investigate the differential capacitance profiles and interfacial structures of MoS₂ electrodes in ILs.

In this work, we studied the energy storage mechanisms of MoS₂-based supercapacitors in imidazolium ionic liquid ([BMI⁺][PF₆⁻]) using MD simulation. [BMI⁺][PF₆⁻] ionic liquid is chosen because of its high electrical stability window (~4.6 V), superior thermal stability (up to 400 °C), and suitable electrical conductivity (~2–5 mS/cm), which has been widely studied as the electrolyte of supercapacitors. The interfacial EDL structure is comprehensively investigated, including number density, charge density, electric potential distributions, electric field distributions, and separation coefficient, which is further explained in terms of ion-electrode interaction energy. Especially, a bell-shaped differential capacitance profile is recognized, which can be well described by the differential charge density curve in the Helmholtz layer. Moreover, the diffusion coefficient properties within different EDL regions are evaluated, revealing that ion dynamic properties are significantly suppressed by ~50% in the Helmholtz region.

2. Simulation Methods

The molecular dynamics simulation model is composed of positively charged MoS₂ (right), negatively charged MoS₂ (left), and [BMI⁺][PF₆⁻] electrolytes in the middle region. The dimensions of MoS₂ electrodes along the X and Y directions were ~46.9646 Å and ~44.66 Å, respectively. In this work, partial charges of 0 μC cm⁻², ±4 μC cm⁻², ±8 μC cm⁻², ±12 μC cm⁻², and ±16 μC cm⁻² were placed on the MoS₂ electrode to mimic the applied electric potentials. Our previous studies evaluated the reliability of the constant charge method for MD simulations in predicting the EDL structure and charging dynamics, suggesting that the two methods yield essentially similar interfacial structures and capacitance performance in the equilibrium state [26,27]. However, the constant charge method has much less computational cost (by one order of magnitude) compared to the constant potential method. As a result, the constant charge method has been widely used for studying the interfacial EDL structures and capacitive behaviors in the equilibrium state [16,20,28,29].

The force fields for the MoS₂ electrode are taken from the simulation work of V. Varshney [30], including charges, bond, angle, and non-bond parameters. A four-site coarse-grained model of [BMI⁺][PF₆⁻] ionic liquid is employed in this work. The force fields for the coarse-grained model of [BMI⁺][PF₆⁻] electrolyte are taken from the simulation work of C. Merlet [31]. The polarizations of electrolytes are neglected in the work. Different from the kinetic properties of ionic liquids (e.g., charging dynamics) that require a polarizable force field [19,22], the majority of electrolyte structural properties can be captured accurately by the non-polarizable force fields [32]. In these models, the interactions between atomic sites (or ions) were expressed as a sum of coulombic and Lennard–Jones 12-6 potentials. The cross parameters are obtained using the Lorentz-Berthelot mixing rules.

MD simulations were performed using the LAMMPS software (version: 16Feb16). The Nose-Hoover thermostat with a 100 fs damping parameter was used to control the temperature of the electrolyte (~400 K). A cutoff of ~10.0 Å was used for the calculation of van der Waals (VDW) and coulombic interactions in real space. The long-range coulombic interactions were treated by the particle–particle particle–mesh algorithm with a root mean square accuracy of 10⁻⁵. The Velocity Verlet algorithm was employed to integrate Newton’s equation with a time step of 1 fs. Periodic boundary conditions are applied along the X and Y direction, and the vacuum layer is set along the Z direction to avoid periodic interactions. During the MD simulation, the electrode atoms of MoS₂ are kept fixed (zero velocity) with additional forces. The simulation of ionic liquids is initially performed in the isobaric−isothermal ensemble (NPT) to obtain equilibrium property, which is then placed into the middle region of two MoS₂ electrodes. The box size along the Z direction (130 Å) is adjusted to yield bulk density in the middle of the simulation box. In the NVT ensemble, the MD system is first heated at 1000 K for 2 ns and then annealed to 400 K over a period of 1 ns, followed by another 7 ns to reach equilibrium. After that, a 20 ns production is performed for data analysis.

3. Results and Discussions

MD simulation model of MoS₂ electrode in [BMI⁺][PF₆⁻] electrolyte is shown in Figure 1, which is used to investigate the interfacial behaviors, dynamics, and charge storage performance. The constant charge method is applied to mimic the electric potential applied on the MoS₂ electrode. Similar MD models and methods can be widely observed in previous theoretical studies [33,34].

Firstly, to study the interfacial structures of MoS₂ electrodes, the number density of [BMI⁺][PF₆⁻] electrolyte at different charge densities is characterized, defined as:

(1)${\rho }_n\left(z\right)=\frac{1}{XY\Delta z}{{\displaystyle \sum }}_i\delta \left(z-z_i\right)$

where ${{\displaystyle \sum }}_i\delta \left(z-z_i\right)$ present the number density of i ion at a slab of z − z_i.

Figure 2a,b show the interfacial populations of BMI⁺ cations and PF₆⁻ anions nearby the MoS₂ electrodes. Based on the number density distribution, we define the layers of EDL structures as follows: (1) inner Helmholtz layer located within 5.0 Å from the electrode surface, (2) outer Helmholtz layer located between 5.0 and 10 Å from the electrode surface, and (3) diffusion layer located > 10 Å from the electrode surface.

On the uncharged MoS₂ electrode, [BMI⁺][PF₆⁻] electrolyte shows a feature of multiple alternating layers of BMI⁺ ions and PF₆⁻ ions at the electrode-ILs interface. This can be explained in terms of strong electrode-ILs van der Waals interactions and cation–anion correlations in solvent-free ILs systems. On the uncharged MoS₂ electrode, the inner Helmholtz layer is mainly composed of BMI⁺ ion and PF₆⁻ ion.

Upon applying negative charges, more BMI⁺ ions are attracted into the inner Helmholtz layer; meanwhile, desorption of PF₆⁻ ions occurs, which are moved away into the outer Helmholtz layer. As a result, the inner Helmholtz layer is mainly composed of BMI⁺ ions, which is followed by a second layer of PF₆⁻ ion in the outer Helmholtz layer. Such a multilayer structure propagates several nanometers away from the electrode, particularly at high electrode charge density.

With increasing the amount of charge density to −16 μC/cm², the magnitude of the first peak of BMI⁺ cation in the inner Helmholtz layer is significantly enhanced by 2.7 times (from 0.01736 #/ų to 0.04672 #/ų). Its packing position of the first peak is also shifted closer to the MoS₂ electrode (from 3.7 Å to 3.4 Å). Meanwhile, the magnitude of the first peak of PF₆⁻ anion is significantly reduced by 52.5% (from 0.02316 #/ų to 0.01098 #/ų), with the charge density increasing to −4 μC/cm². Upon further increasing negative charges, the position of the first peak of the PF₆⁻ anion is moved further away from the MoS₂ electrode (from 4.0 Å to 7.0 Å) into the outer Helmholtz layer. However, further increasing the negative charges, the interfacial anion population in the outer Helmholtz layer is slightly improved (from 0.009 #/ų to 0.014 #/ų). This can be attributed to the intense cation–anion interactions, which are also reported in previous literature [35].

Figure 2c,d display the interfacial populations of BMI⁺ cations and PF₆⁻ anions nearby the positively charged MoS₂ electrodes. Upon increasing the positive charges, the density of the first peak of PF₆⁻ anion in the inner Helmholtz layer is obviously enlarged by 2.75 folds (from 0.0257 #/ų to 0.0708 #/ų) with the location at 3.6 Å to MoS₂ electrode. As expected, the intensity of BMI⁺ cation is decreased (from 0.0175 #/ų to 0.0072 #/ų at 4 μC/cm²), where its packing position is shifted to 6.6 Å into the outer Helmholtz layer, and the number density is improved to 0.0115 #/ų at 16 μC/cm². This section may be divided into subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Representative MD snapshots are shown in Figure 3 to reveal the electrolyte layering structure at different charge densities. Upon increasing negative charges, the interfacial population of cation within the inner Helmholtz region is enhanced; meanwhile, the anion population is suppressed. Especially, no anions can be observed at −16 μC/cm² due to the strong electrostatic repulsion interactions. Similar results can be observed on the positive polarization. Differently, cations can still be observed even at 16 μC/cm², which can be ascribed to the specific adsorption of BMI⁺ on the MoS₂ surface.

To describe the ion–MoS₂ and cation–anion interactions, the separation coefficient of ion within the Helmholtz layer is calculated, which is defined as:

(2)$\theta =\frac{{{\displaystyle \int }}_{z_0}^{z_1}{\rho }_{\mathrm{counterion}}\left(z\right)\mathit{dz}}{{{\displaystyle \int }}_{z_0}^{z_1}{\rho }_{\mathrm{counterion}}\left(z\right)+{\rho }_{\mathrm{coion}}\left(z\right)\mathit{dz}}\times 100\%$

Based on the number density results, the span of the inner Helmholtz region is about 5 Å on the MoS₂ electrode, which is chosen for the separation coefficient calculation (Z₁ = 5 Å). As shown in Figure 4, upon the negative polarization, the separation coefficient of counterion (i.e., BMI⁺) is prominently improved from 62.5% at −4 μC/cm² to 89.5% at −8 μC/cm², 99.9% at −12 μC/cm², and 100% at −16 μC/cm², which agrees well with the number density results. Similar results can be observed on the positive MoS₂ electrode, which is increased to 99.5% at 16 μC/cm².

The interaction energies between ions in the interfacial EDL regions and MoS₂ electrode are evaluated, which helps to explain the ion population in the Helmholtz layer to some extent. As shown in Figure 5, the interaction energy value is around a constant value during the simulation time from 10 ns to 30 ns, indicating that ions nearby the MoS₂ electrode are in an equilibrium state. For the negatively charged MoS₂ electrode, the interaction energy of MoS₂-BMI⁺ is approx. −140 kcal/mol, which is much higher than that of PF₆⁻ anion (−10 kcal/mol). On the contrary, nearby the positively charged MoS₂ electrode, the interaction energy of MoS₂-BMI⁺ (−7 kcal/mol) is much lower than that of PF₆⁻ anion (−210 kcal/mol). These results can account for the formation mechanisms of the inner and outer Helmholtz layer on the MoS₂ electrode.

In MD simulations, the electric potential distributions can be calculated by integrating the Poisson equations:

(3)${\nabla }_z\left[{\epsilon }_0\left({\nabla }_z\phi \left(z\right)\right)\right]=-\sigma \left(z\right)$

where ${\mathsf{\epsilon }}_0$, $\sigma \left(z\right)$, and $\phi \left(z\right)$ show the permittivity of vacuum, the charge density at z, and the electric potential at z, respectively.

The electric potential distribution across the simulation cell can be obtained:

(4)$\phi (z)=\frac{\alpha }{{\epsilon }_0}z-\frac{1}{{\epsilon }_0}{\displaystyle {\int }_0^z(z-z^{\prime })\sigma (z^{\prime })d{z}^{\prime }}$

Considering the complete screening of interfacial electric fields in the bulk region, the electrode potential U_electrode was calculated, defined as the potential drop within the electrode/electrolyte interfacial district relative to the potential of zero charge (PZC):

(5)$U_{\mathrm{electrode}}=U_{\mathrm{EDL}}-\mathit{PZC}={\phi }_{\mathrm{electrode}}-{\phi }_{\mathrm{bulk}}-\mathit{PZC}$

where U_EDL is defined as the potential drop across the EDL-dominant area; PZC is the potential of zero charge.

Figure 6a presents the electric potential distributions nearby the negative MoS₂ electrode. A multilayer structure feature is also recognized, especially at high charge density, which is in suitable accordance with the number density distributions. Especially, the peak position of the electric potential curve can correspond to the location of BMI⁺ ion populations. For example, the electric potential profile at −16 μC/cm² presents a first peak at 3.6 Å, a second peak at 10.5 Å, and a third peak at 17 Å, which is in suitable agreement with BMI⁺ ion number density. Meanwhile, the electric potential profile at −16 μC/cm² presents a first valley at 7.2 Å and a second valley at 13.8 Å, which can be indexed to the peak position of PF₆⁻ ion number density. Interestingly, we would like to highlight that the correlations between the peak/valley position of the electric potential curve and cation/anion number density profiles have never been reported before. As shown in Figure 6a, the uncharged MoS₂ electrode shows a PZC value of 0.016 V, which is probably due to the specific adsorption of BMI⁺ ions. With the increment of electrode charge density, the electric potential drop across the EDL region increases from 0.88 V at −4 μC/cm² to 4.22 V at −16 μC/cm².

The interfacial electric field distribution is also presented to understand the screening of electrode charges by electrolyte structure (Figure 6b). It is found that the strength of the electric field is rapidly suppressed in the inner and outer Helmholtz layers. Especially, a more obvious multilayer structure is recognized. It shows a first peak at 8.2 Å, a second peak at 15.1 Å, and a third peak at 21.8 Å; meanwhile, it shows the first valley at 5.3 Å, the second valley at 12.2 Å, and the third valley at 18.8 Å. These features basically agree well with the electrolyte charge density distributions (Figure 6c). The position of the peak corresponds to the position of negative charges, while the valley location can be indexed to the position of positive charges. These results help to understand the screening of electrode charges by interfacial EDL structures. Similar results can be observed on the positive MoS₂ electrodes (Figure 6d–f).

To qualitatively study the response of EDL structure upon the applied electric fields, we calculate the differential capacitance curves of MoS₂ electrodes in [BMI⁺][PF₆⁻] electrolyte. The differential capacitance can be calculated by:

(6)$C_D=\frac{d\alpha }{d{U}_{\mathrm{EDL}}}$

where α is the charge density of the MoS₂ electrode and U_EDL is the electric potential drop across the EDL region.

As shown in Figure 7a, the charge density increases monotonically with enlarging the electric potentials, which is in an almost linear relationship. Especially, the MoS₂ electrode shows bell-shaped differential capacitance profiles with respect to the electric potentials (Figure 7b), which is different from the U-shaped curve from conventional Gouy–Chapman theory. A maximum capacitance of 4.86 μF/cm² is observed at 0 V (with respect to PZC value), which is decreased to 1.23 and 2.94 μF/cm² upon the negative and positive polarization, respectively.

The bell-shaped differential capacitance curves are closely correlated with the high ion density at the interface, as observed in Figure 2, especially within the inner Helmholtz layer. The high population of cations and anions nearby the uncharged MoS₂ electrode induce strong excluded volume steric interactions that are closely related to the compressibility. This will induce a high free-energy penalty for ion transporting from the bulk region to the interface to form EDL structures, which results in a much slower extent of counterion accumulation near the MoS₂ electrode surface with the increment of electric potential, leading to reduced capacitive values. Kornyshev et al. also demonstrated that the volume steric interactions and electrolyte compressibility within the EDL region play a major role in determining differential capacitance shape [36]. The as-obtained capacitance values are also close to previous theoretical and experimental studies [37,38]. For example, J. Vatamanu suggested that the differential capacitance values of ILs (e.g., [pyr₁₃⁺][FSI⁻] and [pyr₁₃⁺][TFSI⁻]) are around 2~ 5 μF/cm² at different electric potentials [39].

According to the electric double-layer theory, the specific capacitance of a single electrode could be calculated as:

(7)$C_{\mathrm{electrode}}=\frac{\alpha }{U_{\mathrm{electrode}}}$

where α is the surface charge density of electrodes.

Figure 7c shows the electrode capacitance at different surface charge densities. At the relatively low charge density, the capacitance value at the negative electrode is higher than that of the positive counterpart. For example, the capacitance value of C⁻ (~4.63 μF/cm²) is much larger than that of C⁺ (~3.88 μF/cm²) at 4 μC/cm². However, such a trend decreases with the increment of charge density, and an opposite phenomenon is observed at 16 μC/cm², where the capacitance value of C⁻ (~3.81 μF/cm²) is slightly lower than that of C⁺ (~4.06 μF/cm²). Similar results are also reported on the graphene electrodes. For example, C. Merlet reported that the electrode capacitance of [BMI⁺][PF₆⁻] is 4.8 and 3.8 μF/cm² on the negative and positive graphene electrodes, respectively [31].

To establish the relationship between interfacial EDL structures and capacitive behaviors, the differential charge density is introduced and calculated by:

(8)$C_D=\frac{d{\alpha }_e}{d{U}_{\mathrm{EDL}}}$

where α_e is the cumulated charge density of electrolyte in the EDL region, and U_EDL is the electric potential drop across the EDL region. Cumulated electrolyte charge densities in different EDL regions at the distance of 5, 7.5, 10, 12.5, 15, 17.5, and 20 Å to the MoS₂ electrode are calculated. As shown in Figure 8, differential charge density profiles within different EDL regions show completely different shapes. Specifically, the differential charge density profiles at the distance of 5, 10, 12.5, and 17.5 Å to the MoS₂ electrode exhibit bell-shaped curves; the differential charge density profiles at the distance of 7.5, 15, and 20 Å to the MoS₂ electrode present the asymmetric U-shaped curves. In aqueous electrolytes, the differential charge density curve can well reproduce that of differential capacitance profiles in the 10 Å-thickness EDL region. However, in ILs systems, the shape of the differential charge density curve in the 5 Å thickness EDL region is very close to the belled-shape differential capacitance profiles. This reveals that the capacitive behaviors of [BMI⁺][PF₆⁻] electrolyte are dominantly determined by the first peak of counterion in the Helmholtz layer (5 Å to the MoS₂ electrode surface), while the electrolyte structure in the outer Helmholtz layer and diffusion layer plays a minor role. This finding also indicates that the effective thickness of the EDL region is approx. 5 Å.

In addition to the interfacial layering structures, the dynamic properties of [BMI⁺][PF₆⁻] electrolyte in the EDL region with different distances to the MoS₂ electrode surface are studied, which helps to understand the charging mechanisms of ILs-based supercapacitors. The self-diffusion coefficients D are calculated based on the mean square displacement method:

(9)$D=\frac{1}{2n}\underset{t\to \infty }{\lim }\frac{{⟨r\left(t\right)-r\left(0\right)⟩}^2}{t}$

where r(t) is the position of the particle mass center at the time t, and n is the dimensionality of the system. Ion diffusion coefficients in the EDL region located at 0~5 Å, 5~10 Å, and 10~15 Å from MoS₂ electrodes are calculated, which corresponds to the inner Helmholtz layer, outer Helmholtz layer, and diffusion layer, respectively. The dynamic properties in the bulk electrolyte far away from the interfacial region are also evaluated for comparison.

As shown in Figure 9, the diffusion coefficient of BMI⁺ cation in the bulk region is 3.33 × 10⁻⁹ m² s⁻¹. The ion dynamics are significantly reduced in the interfacial region of MoS₂ electrodes. Specifically, the diffusion coefficients of BMI⁺ cation in the EDL region of 0~5 Å and 5~10 Å distance to the MoS₂ electrode surface are 1.32 × 10⁻⁹ m² s⁻¹ and 1.53 × 10⁻⁹ m² s⁻¹, respectively, which is reduced by ~60.0% and ~54.5%, respectively, compared to the bulk state. Moreover, the diffusion coefficient increases significantly with moving away from the MoS₂ electrode, which is calculated to be 3.26 × 10⁻⁹ m² s⁻¹ in the EDL region of 10~15 Å from the electrode surface. Especially, the corresponding diffusion properties are close to the bulk properties with a distance over 10 Å. The observed reduced kinetics are mainly attributed to the strong ion-electrode and ion-ion interactions in the inner and outer Helmholtz layers. Similar results are observed in the positively charged systems. The diffusion coefficients of PF₆⁻ anions are calculated to be 1.3 × 10⁻⁹ m² s⁻¹, 2.19 × 10⁻⁹ m² s⁻¹, and 2.77 × 10⁻⁹ m² s⁻¹ in the interfacial region of 0~5 Å, 5~10 Å, and 10~15 Å from the positive MoS₂ electrode, respectively. In the inner and outer Helmholtz layers, the dynamics are reduced by ~53.7% and ~22.3%, respectively, compared to the bulk dynamics (2.82 × 10⁻⁹ m² s⁻¹).

The MD simulation trajectories of BMI⁺ ion and PF₆⁻ ion are also recorded to visualize the cation–anion correlations and dynamic properties. As shown in Figure 10a, the trajectories of BMI⁺ ion and PF₆⁻ ion are very close to each other during the simulation. This indicates the strong correlations between BMI⁺ ion and PF₆⁻ ion, which is mainly in the form of electrostatic interactions. The intense ion-ion correlations can be widely observed in the solvent-free ILs systems, which become much weaker in the organic and aqueous electrolytes.

Radical distribution function g(r) is calculated to quantify the correlations between BMI⁺ ion and PF₆⁻ ion. It is found that g(r) between BMI⁺ ion and PF₆⁻ ion shows an obvious peak at 4.85 Å (peak magnitude of 5.97). The corresponding coordination number is calculated to be 1.3, revealing the existence of cation-anions pairs.

4. Conclusions

In summary, the atomistic-level energy storage mechanisms of MoS₂ electrode in imidazolium ionic liquid ([BMI⁺][PF₆⁻]) are explored using MD simulation. The EDL structures of MoS₂ electrodes in [BMI⁺][PF₆⁻] electrolytes are comprehensively studied in terms of number density, charge density, electric potential distribution, electric field distribution, and separation coefficient. The EDL structure can be explained in terms of ion-electrode interaction energy. A bell-shaped differential capacitance profile is recognized, which can be reproduced by the differential charge density curve in the inner Helmholtz layer. This indicates that the capacitive behaviors of the MoS₂ electrode are mainly determined by the counterion population/structure within the 5.0 Å EDL region from the electrode surface. Moreover, the diffusion coefficient properties within different EDL regions are evaluated, revealing that ion dynamic properties are significantly suppressed by ~50% in the Helmholtz region, which can be recovered to a bulk state with a location 10 Å away from the electrode. The as-obtained insights can provide useful guidance on improving the energy density of MoS₂ supercapacitors. Future studies should be paid to investigating the charging dynamics and capacitive performance of nanometer MoS₂ pores in ionic liquids, which help to boost the practical applications of MoS₂-based supercapacitors.

Author Contributions

Conceptualization, C.X., Z.X. and H.Y.; methodology, C.X.; investigation, C.X., Z.X., Y.W., J.Y., H.C. and Q.L.; formal analysis, C.X., Z.X. and J.Y.; visualization, C.X.; validation, H.C., Q.L. and G.C.; supervision, H.Y.; resources, H.Y.; writing—original draft, C.X.; writing—review and editing, Y.W., J.Y., H.C., Q.L., G.C. and H.Y.; funding acquisition, H.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. MD simulation model of [BMI+][PF6−] electrolyte on the MoS2 electrode. Blue and orange colors present the cation and anion, respectively. Green and yellow colors denote Mo atom and S atom, respectively.

Figure 2. Number density of (a) BMI+ cation and (b) PF6− anion nearby the negatively charged MoS2 electrode. Number density of (c) BMI+ cation and (d) PF6− anion nearby the positively charged MoS2 electrode.

Figure 3. Representative MD simulation snapshot of [BMI+][PF6−] electrolyte on the MoS2 electrodes. Green and yellow colors denote the electrode of the Mo atom and S atom, respectively. Purple and lemon-yellow colors present the electrolyte of cation and anion, respectively.

Figure 4. Separation coefficient at different surface charge densities.

Figure 5. Interaction energy between (a) BMI+, (b) PF6− and negatively charged MoS2 electrode. Interaction energy between (c) BMI+, (d) PF6− and positively charged MoS2 electrode.

Figure 6. (a) Electric potential, (b) electric field, and (c) electrolyte charge density profile nearby negatively charged MoS2 electrode. (d) Electric potential, (e) electric field, and (f) electrolyte charge density profile nearby positively charged MoS2 electrode.

Figure 7. (a) Charge density as a function of electric potential. (b) Differential capacitance curves. (c) Electrode capacitance as a function of charge density.

Figure 8. Differential electrolyte charge density profile as a function of electric potential.

Figure 9. (a) Diffusion coefficients of BMI+ ions in different EDL regions on the negative MoS2 electrode. (b) Diffusion coefficients of PF6− ions in different EDL regions on the positive MoS2 electrode.

Figure 10. (a) MD simulation trajectory of BMI+ ion (red color) and PF6− ion (blue color). (b) Radical distribution function between BMI+ ion and PF6− ion.

References

1. Simon, P.; Gogotsi, Y. Perspectives for electrochemical capacitors and related devices. Nat. Mater.; 2020; 19, pp. 1151-1163. [DOI: https://dx.doi.org/10.1038/s41563-020-0747-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32747700]

2. Wang, J.; Malgras, V.; Sugahara, Y.; Yamauchi, Y. Electrochemical energy storage performance of 2D nanoarchitectured hybrid materials. Nat. Commun.; 2021; 12, 3563. [DOI: https://dx.doi.org/10.1038/s41467-021-23819-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34117228]

3. Xiao, J.; Han, J.; Zhang, C.; Ling, G.; Kang, F.; Yang, Q.-H. Dimensionality, Function and Performance of Carbon Materials in Energy Storage Devices. Adv. Energy Mater.; 2022; 12, 2100775. [DOI: https://dx.doi.org/10.1002/aenm.202100775]

4. Pomerantseva, E.; Bonaccorso, F.; Feng, X.; Cui, Y.; Gogotsi, Y. Energy storage: The future enabled by nanomaterials. Science; 2019; 366, eaan8285. [DOI: https://dx.doi.org/10.1126/science.aan8285]

5. Liu, L.; Taberna, P.-L.; Dunn, B.; Simon, P. Future Directions for Electrochemical Capacitors. ACS Energy Lett.; 2021; 6, pp. 4311-4316. [DOI: https://dx.doi.org/10.1021/acsenergylett.1c01981]

6. Nasrin, K.; Sudharshan, V.; Subramani, K.; Sathish, M. Insights into 2D/2D MXene Heterostructures for Improved Synergy in Structure toward Next-Generation Supercapacitors: A Review. Adv. Funct. Mater.; 2022; 32, 2110267. [DOI: https://dx.doi.org/10.1002/adfm.202110267]

7. Voiry, D.; Fullon, R.; Yang, J.; de Carvalho Castro e Silva, C.; Kappera, R.; Bozkurt, I.; Kaplan, D.; Lagos, M.J.; Batson, P.E.; Gupta, G. et al. The role of electronic coupling between substrate and 2D MoS₂ nanosheets in electrocatalytic production of hydrogen. Nat. Mater.; 2016; 15, pp. 1003-1009. [DOI: https://dx.doi.org/10.1038/nmat4660]

8. Jiao, Y.; Hafez, A.M.; Cao, D.; Mukhopadhyay, A.; Ma, Y.; Zhu, H. Metallic MoS₂ for High Performance Energy Storage and Energy Conversion. Small; 2018; 14, 1800640. [DOI: https://dx.doi.org/10.1002/smll.201800640]

9. Geng, X.; Jiao, Y.; Han, Y.; Mukhopadhyay, A.; Yang, L.; Zhu, H. Freestanding Metallic 1T MoS₂ with Dual Ion Diffusion Paths as High Rate Anode for Sodium-Ion Batteries. Adv. Funct. Mater.; 2017; 27, 1702998. [DOI: https://dx.doi.org/10.1002/adfm.201702998]

10. Acerce, M.; Voiry, D.; Chhowalla, M. Metallic 1T phase MoS₂ nanosheets as supercapacitor electrode materials. Nat. Nanotechnol.; 2015; 10, pp. 313-318. [DOI: https://dx.doi.org/10.1038/nnano.2015.40]

11. Feng, N.; Meng, R.; Zu, L.; Feng, Y.; Peng, C.; Huang, J.; Liu, G.; Chen, B.; Yang, J. A polymer-direct-intercalation strategy for MoS₂/carbon-derived heteroaerogels with ultrahigh pseudocapacitance. Nat. Commun.; 2019; 10, 1372. [DOI: https://dx.doi.org/10.1038/s41467-019-09384-7]

12. Bo, Z.; Cheng, X.; Yang, H.; Guo, X.; Yan, J.; Cen, K.; Han, Z.; Dai, L. Ultrathick MoS₂ Films with Exceptionally High Volumetric Capacitance. Adv. Energy Mater.; 2022; 12, 2103394. [DOI: https://dx.doi.org/10.1002/aenm.202103394]

13. Ye, J.; Simon, P.; Zhu, Y. Designing ionic channels in novel carbons for electrochemical energy storage. Natl. Sci. Rev.; 2020; 7, pp. 191-201. [DOI: https://dx.doi.org/10.1093/nsr/nwz140]

14. Jeanmairet, G.; Rotenberg, B.; Salanne, M. Microscopic Simulations of Electrochemical Double-Layer Capacitors. Chem. Rev.; 2022; 122, pp. 10860-10898. [DOI: https://dx.doi.org/10.1021/acs.chemrev.1c00925]

15. Yang, H.; Bo, Z.; Yang, J.; Kong, J.; Chen, X.; Yan, J.; Cen, K. Substrate Effects in Graphene-Based Electric Double-Layer Capacitors: The Pivotal Interplays between Ions and Solvents. ChemElectroChem; 2017; 4, pp. 2966-2974. [DOI: https://dx.doi.org/10.1002/celc.201700733]

16. Yang, H.; Bo, Z.; Yan, J.; Cen, K. Influence of wettability on the electrolyte electrosorption within graphene-like nonconfined and confined space. Int. J. Heat Mass Transf.; 2019; 133, pp. 416-425. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.12.134]

17. Yang, H.; Yang, J.; Li, C.; Huang, Z.; Bendavid, A.; Yan, J.; Cen, K.; Han, Z.; Bo, Z. Gel polymer dominated ion charging mechanisms within graphene nanochannels. J. Power Sources; 2022; 541, 231684. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2022.231684]

18. Seebeck, J.; Schiffels, P.; Schweizer, S.; Hill, J.-R.; Meißner, R.H. Electrical Double Layer Capacitance of Curved Graphite Electrodes. J. Phys. Chem. C; 2020; 124, pp. 5515-5521. [DOI: https://dx.doi.org/10.1021/acs.jpcc.9b10428]

19. Tu, Y.-J.; McDaniel, J.G. Structure–Capacitance Relationships of Graphene/Ionic Liquid Electrolyte Double Layers. J. Phys. Chem. C; 2021; 125, pp. 20204-20218. [DOI: https://dx.doi.org/10.1021/acs.jpcc.1c05698]

20. Chaban, V.V.; Andreeva, N.A.; Fileti, E.E. Graphene/ionic liquid ultracapacitors: Does ionic size correlate with energy storage performance?. New J. Chem.; 2018; 42, pp. 18409-18417. [DOI: https://dx.doi.org/10.1039/C8NJ04399J]

21. Bi, S.; Wang, R.; Liu, S.; Yan, J.; Mao, B.; Kornyshev, A.A.; Feng, G. Minimizing the electrosorption of water from humid ionic liquids on electrodes. Nat. Commun.; 2018; 9, 5222. [DOI: https://dx.doi.org/10.1038/s41467-018-07674-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30514881]

22. Tu, Y.-J.; Delmerico, S.; McDaniel, J.G. Inner Layer Capacitance of Organic Electrolytes from Constant Voltage Molecular Dynamics. J. Phys. Chem. C; 2020; 124, pp. 2907-2922. [DOI: https://dx.doi.org/10.1021/acs.jpcc.0c00299]

23. Zhan, C.; Lian, C.; Zhang, Y.; Thompson, M.W.; Xie, Y.; Wu, J.; Kent, P.R.C.; Cummings, P.T.; Jiang, D.-E.; Wesolowski, D.J. Computational Insights into Materials and Interfaces for Capacitive Energy Storage. Adv. Sci.; 2017; 4, 1700059. [DOI: https://dx.doi.org/10.1002/advs.201700059] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28725531]

24. Sun, C.; Zhou, R.; Zhao, Z.; Bai, B. Nanoconfined Fluids: What Can We Expect from Them?. J. Phys. Chem. Lett.; 2020; 11, pp. 4678-4692. [DOI: https://dx.doi.org/10.1021/acs.jpclett.0c00591]

25. Xu, K.; Shao, H.; Lin, Z.; Merlet, C.; Feng, G.; Zhu, J.; Simon, P. Computational Insights into Charge Storage Mechanisms of Supercapacitors. Energy Environ. Mater.; 2020; 3, pp. 235-246. [DOI: https://dx.doi.org/10.1002/eem2.12124]

26. Yang, H.; Yang, J.; Bo, Z.; Chen, X.; Shuai, X.; Kong, J.; Yan, J.; Cen, K. Kinetic-Dominated Charging Mechanism within Representative Aqueous Electrolyte-based Electric Double-Layer Capacitors. J. Phys. Chem. Lett.; 2017; 8, pp. 3703-3710. [DOI: https://dx.doi.org/10.1021/acs.jpclett.7b01525]

27. Yang, J.; Bo, Z.; Yang, H.; Qi, H.; Kong, J.; Yan, J.; Cen, K. Reliability of Constant Charge Method for Molecular Dynamics Simulations on EDLCs in Nanometer and Sub-Nanometer Spaces. ChemElectroChem; 2017; 4, pp. 2486-2493. [DOI: https://dx.doi.org/10.1002/celc.201700447]

28. Begić, S.; Jónsson, E.; Chen, F.; Forsyth, M. Molecular dynamics simulations of pyrrolidinium and imidazolium ionic liquids at graphene interfaces. Phys. Chem. Chem. Phys.; 2017; 19, pp. 30010-30020. [DOI: https://dx.doi.org/10.1039/C7CP03389C]

29. Wu, M.; Li, W.; Li, S.; Feng, G. Capacitive performance of amino acid ionic liquid electrolyte-based supercapacitors by molecular dynamics simulation. RSC Adv.; 2017; 7, pp. 28945-28950. [DOI: https://dx.doi.org/10.1039/C7RA00443E]

30. Varshney, V.; Patnaik, S.S.; Muratore, C.; Roy, A.K.; Voevodin, A.A.; Farmer, B.L. MD simulations of molybdenum disulphide (MoS₂): Force-field parameterization and thermal transport behavior. Comput. Mater. Sci.; 2010; 48, pp. 101-108. [DOI: https://dx.doi.org/10.1016/j.commatsci.2009.12.009]

31. Merlet, C.; Salanne, M.; Rotenberg, B.; Madden, P.A. Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations. J. Phys. Chem. C; 2011; 115, pp. 16613-16618. [DOI: https://dx.doi.org/10.1021/jp205461g]

32. Vatamanu, J.; Hu, Z.; Bedrov, D.; Perez, C.; Gogotsi, Y. Increasing Energy Storage in Electrochemical Capacitors with Ionic Liquid Electrolytes and Nanostructured Carbon Electrodes. J. Phys. Chem. Lett.; 2013; 4, pp. 2829-2837. [DOI: https://dx.doi.org/10.1021/jz401472c]

33. Bo, Z.; Xu, C.; Huang, Z.; Chen, P.; Yan, G.; Yang, H.; Yan, J.; Cen, K.; Ostrikov, K.K. Photo-electric capacitive deionization enabled by solar-driven nano-ionics on the edges of plasma-made vertical graphenes. Chem. Eng. J.; 2021; 422, 130156. [DOI: https://dx.doi.org/10.1016/j.cej.2021.130156]

34. Xu, C.; Zhu, J.; Li, D.; Xu, Q.; Chen, G.; Yang, H. Unveiling the Effects of Solvent Polarity within Graphene Based Electric Double-Layer Capacitors. Energies; 2022; 15, 9487. [DOI: https://dx.doi.org/10.3390/en15249487]

35. Silvester, D.S.; Jamil, R.; Doblinger, S.; Zhang, Y.; Atkin, R.; Li, H. Electrical Double Layer Structure in Ionic Liquids and Its Importance for Supercapacitor, Battery, Sensing, and Lubrication Applications. J. Phys. Chem. C; 2021; 125, pp. 13707-13720. [DOI: https://dx.doi.org/10.1021/acs.jpcc.1c03253]

36. Kornyshev, A.A. Double-Layer in Ionic Liquids:  Paradigm Change?. J. Phys. Chem. B; 2007; 111, pp. 5545-5557. [DOI: https://dx.doi.org/10.1021/jp067857o]

37. Lamperski, S.; Outhwaite, C.W.; Bhuiyan, L.B. The Electric Double-Layer Differential Capacitance at and near Zero Surface Charge for a Restricted Primitive Model Electrolyte. J. Phys. Chem. B; 2009; 113, pp. 8925-8929. [DOI: https://dx.doi.org/10.1021/jp900037h]

38. Su, Y.-Z.; Fu, Y.-C.; Yan, J.-W.; Chen, Z.-B.; Mao, B.-W. Double Layer of Au(100)/Ionic Liquid Interface and Its Stability in Imidazolium-Based Ionic Liquids. Angew. Chem. Int. Ed.; 2009; 48, pp. 5148-5151. [DOI: https://dx.doi.org/10.1002/anie.200900300]

39. Vatamanu, J.; Borodin, O.; Smith, G.D. Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-N-propylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes. J. Phys. Chem. B; 2011; 115, pp. 3073-3084. [DOI: https://dx.doi.org/10.1021/jp2001207]