Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode Knowledge of how the molecular structures of ionic liquids (ILs) affect their properties at electrified interfaces is key towards the rational design of ILs for electric applications. Polarizable molecular dynamics simulations were performed to investigate the structural, electrical and dynamic properties of electric double layers (EDLs) formed by imidazolium dicyanamide ([ImX1][DCA]) at the interface with molybdenum disulfide electrode. The effect of side chain of imidazolium on the properties of EDLs was analyzed by using 1-ethyl-3-methylimidazolium ([Im21]), 1-octyl-3-methylimidazolium ([Im81]), 1-benzyl-3-methylimidazolium ([ImB1]) and 1-(2-hydroxyethyl)-3-methylimidazolium ([ImO1]) as cations. Using [Im21] as reference, we find that the introduction of octyl or benzyl groups significantly alters the interfacial structures near the cathode because of the reorientation of cations. For [Im81], the positive charge on the cathode induces pronounced polar and non-polar domain separation. In contrast, the hydroxyl group has minor effect on the interfacial structures. [ImB1] is shown to deliver slightly larger capacitance than other ILs even though it has larger molecular volume than [Im21]. This is attributed to the limiting factor for capacitance being the strong association between counter-ions, instead of the free space available to ions at the interface. For [Im81], the charging mechanism is mainly the exchange between anions and octyl tails, while for the other ILs, the mechanism is mainly the exchange of counter-ions. Analysis on the charging process shows that the charging speed does not correlate strongly with macroscopic bulk dynamics like viscosity. Instead, it is dominated by local displacement and reorientation of ions.

I. INTRODUCTION

Two-dimensional (2D) nanomaterials such as graphene, hexagonal boron nitride and transition metal dichalcogenide are composed of atomically-thin layers assembled by means of weak, non-covalent forces. Due to electron confinement in two dimensions, 2D materials are able to access unprecedented physical, electronic and chemical properties, therefore be compelling in applications of novel electronic devices and sensors. [START_REF] Geim | The rise of graphene[END_REF][START_REF] Tan | Recent advances in ultrathin two-dimensional nanomaterials[END_REF] In electrical applications such as batteries, electrical double layer capacitors (EDLCs) and fieldeffect transistors (FETs), novel electrolytes are required to function together with 2D nanomaterials. Roomtemperature ionic liquids (ILs) have unique properties, such as wide electrochemical window, high thermal stability, low vapor pressure, making them promising electrolytes. [START_REF] Watanabe | Application of Ionic Liquids to Energy Storage and Conversion Materials and Devices[END_REF][START_REF] Macfarlane | Energy applications of ionic liquids[END_REF][START_REF] Armand | Ionic-liquid materials for the electrochemical challenges of the future[END_REF][START_REF] Fujimoto | Electric-double-layer field-effect transistors with ionic liquids[END_REF][START_REF] Béguin | Carbons and Electrolytes for Advanced Supercapacitors[END_REF] At the electrode-electrolyte interfaces, the ILs form ultra-thin electrical double layers (EDLs). The structure and dynamics of EDLs are among the key factors determining the capacitance and charging speed of EDLC and the gating speed of FET, and thus the performance of the devices. A fundamental understanding of the behavior of EDLs from the molecular level is essential for designing electrolytes with improved performance. In recent decades, numerous research efforts have been reported a) Electronic mail: agilio.padua@ens-lyon.fr on EDLs formed at the IL-electrode interfaces, consisting in experimental, [START_REF] Baldelli | Surface Structure at the Ionic Liquid-Electrified Metal Interface[END_REF] mean-field theory [START_REF] Fedorov | Ionic Liquids at Electrified Interfaces[END_REF] and simulation [START_REF] Vatamanu | Charge storage at the nanoscale: understanding the trends from the molecular scale perspective[END_REF][START_REF] Burt | A review of molecular modelling of electric double layer capacitors[END_REF][START_REF] Merlet | Computer simulations of ionic liquids at electrochemical interfaces[END_REF] studies.

Molecular dynamics (MD) simulation is a powerful tool for investigating the properties of EDLs. It has provided major insights into the structure of the interfacial layers with atomic resolution, helped understand voltagecapacitance relations and the charging mechanisms. [START_REF] Vatamanu | Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-Npropylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes[END_REF][START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF][START_REF] Péan | On the Dynamics of Charging in Nanoporous Carbon-Based Supercapacitors[END_REF][START_REF] Burt | Capacitance of Nanoporous Carbon-Based Supercapacitors Is a Trade-Off between the Concentration and the Separability of the Ions[END_REF][START_REF] Pean | Understanding the different (dis)charging steps of supercapacitors: influence of potential and solvation[END_REF][START_REF] Li | Electric double layer structure and capacitance of imidazolium-based ionic liquids with FSI-and Tf-anions at graphite electrode by molecular dynamic simulations[END_REF] Extensive works have focused on how the properties of EDLs will be affected by the nature of the electrodes, e.g. the curvature, [START_REF] Feng | The importance of ion size and electrode curvature on electrical double layers in ionic liquids[END_REF][START_REF] Feng | Molecular Insights into Carbon Nanotube Supercapacitors: Capacitance Independent of Voltage and Temperature[END_REF] the roughness [START_REF] Vatamanu | On the Influence of Surface Topography on the Electric Double Layer Structure and Differential Capacitance of Graphite/Ionic Liquid Interfaces[END_REF][START_REF] Vatamanu | Increasing Energy Storage in Electrochemical Capacitors with Ionic Liquid Electrolytes and Nanostructured Carbon Electrodes[END_REF] and the presence of nanopores. [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF][START_REF] Péan | On the Dynamics of Charging in Nanoporous Carbon-Based Supercapacitors[END_REF] Notably, it has been shown that nanoporous carbon electrodes could deliver much larger capacitance than planar graphite electrode because the confinement of ions inside the pores suppresses the overscreening effect and allows ions to approach the electrode surface more closely. [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF] The molecular structures of electrolytes is of comparable importance to that of electrodes in determining the properties of EDLs, e.g. imidazolium-and pyridiniumbased ILs are shown to provide larger capacitance but lesser electrochemical stability than ammonium-based ILs. [START_REF] Mousavi | Ionic Liquids as Electrolytes for Electrochemical Double-Layer Capacitors: Structures that Optimize Specific Energy[END_REF] The functionalization of ions offers very broad possibilities to tune the properties of ILs at interfaces. For a planar electrode, it is usually believed that ILs with smaller ion size are likely to deliver larger interfacial capacitance, because of the more compact packing near the electrodes and thinner EDLs, as shown both by experiments [START_REF] Mousavi | Ionic Liquids as Electrolytes for Electrochemical Double-Layer Capacitors: Structures that Optimize Specific Energy[END_REF][START_REF] Lockett | Differential Capacitance of the Electrical Double Layer in Imidazolium-Based Ionic Liquids: Influence of Potential, Cation Size, and Temperature[END_REF] and simulations. [START_REF] Jo | Effects of Alkyl Chain Length on Interfacial Structure and Differential Capacitance in Graphene Supercapacitors: A Molecular Dynamics Simulation Study[END_REF] However, it was recently reported that bulky surface-active an-ions, such as 1,4-bis(2-ethylhexoxy)-1,4-dioxobutane-2sulfonate and dodecyl sulfate, actually deliver larger capacitances than tetrafluoroborate at high temperature, [START_REF] Mao | Self-assembled nanostructures in ionic liquids facilitate charge storage at electrified interfaces[END_REF] and this result was obtained using the same 1-butyl-3methylimidazolium cation. This enhancement was attributed to self-assembly of the surface active anions into lamellar structures that suppress the over-screening effect. This finding emphasizes that the microscopic structures of ions can have a prominent effect on the properties of EDLs, and thus may affect the performance of electric devices dramatically.

In order to rationally design ILs for specific electric applications, it is therefore of particular interest to understand how modifications in the side chains will affect the properties of EDLs. So far, most of the reports focused on the effect of different lengths of alkyl tails, [START_REF] Lockett | Differential Capacitance of the Electrical Double Layer in Imidazolium-Based Ionic Liquids: Influence of Potential, Cation Size, and Temperature[END_REF][START_REF] Jo | Effects of Alkyl Chain Length on Interfacial Structure and Differential Capacitance in Graphene Supercapacitors: A Molecular Dynamics Simulation Study[END_REF][START_REF] Fedorov | Double layer in ionic liquids: The nature of the camel shape of capacitance[END_REF][START_REF] Wang | Molecular Dynamic Simulations of Ionic Liquids at Graphite Surface[END_REF][START_REF] Su | Electric Double Layer of Au(100)/Imidazolium-Based Ionic Liquids Interface: Effect of Cation Size[END_REF][START_REF] Gallegos | Sidechain effects on the capacitive behaviour of ionic liquids in microporous electrodes[END_REF][START_REF] Katakura | Effect of Switching the Length of Alkyl Chains on Electric Double Layer Structure and Differential Capacitance at the Electrode Interface of Quaternary Ammonium-Based Ionic Liquids Studied Using Molecular Dynamics Simulation[END_REF] suggesting that longer side chain results in smaller capacitance and triggers the change of the voltage-capacitance curve from bell-shape to camel-shape. Other than the alkyl tail, the effect of introducing ether groups (increasing the polarity of the side chains) has been considered, showing that the tail containing ether groups is more likely to be parallel to the electrode surface in the first layer than alkyl tails, thus indicating stronger adsorption between ether-functionalized tail on the electrode. [START_REF] Sharma | Structure of Quaternary Ammonium Ionic Liquids at Interfaces: Effects of Cation Tail Modification with Isoelectronic Groups[END_REF] In this work, we investigate the effects of different sizes and chemical functions in side chains upon the local ordering and properties of EDLs MoS 2 is used as electrode here because of its promising application in FETs based on 2D nanomaterials. [START_REF] Perera | Improved Carrier Mobility in Few-Layer MoS2 Field-Effect Transistors with Ionic-Liquid Gating[END_REF][START_REF] Ponomarev | Ambipolar Light-Emitting Transistors on Chemical Vapor Deposited Monolayer MoS2[END_REF][START_REF] Wang | Ionic Liquid Gating of Suspended MoS2 Field Effect Transistor Devices[END_REF][START_REF] Tang | Realizing Wafer-Scale and Low-Voltage Operation MoS2 Transistors via Electrolyte Gating[END_REF] In IL-gated FETs, the electrostatic field exerted by the EDL induces doping on the semiconducting flake MoS 2 electrodes by stabilizing charge carriers with a density up to 10 14 cm -2 . [START_REF] Zhang | Ambipolar MoS2 thin flake transistors[END_REF] A high capacitance at the electrolytic interface is linked to a high carrier density in the semicon-ductor, albeit with a lower carrier mobility. [START_REF] Xia | Carrier Localization on Surfaces of Organic Semiconductors Gated with Electrolytes[END_REF][START_REF] Li | Role of Percolation in the Conductance of Electrolyte-Gated ${\mathrm{SrTiO}}_{3}$[END_REF] These quantities are strongly dependent on the nature of the ILs, which is an argument towards optimizing the electrolyte through judicious choice of the chemical structure of the ions. Recently, the gating effect of amorphous indium-gallium-zinc oxide FET coupled with a series of imidazolium-based ILs was investigated, demonstrating that both charge density and homogeneity of charge distribution at the electrode surface are important to the gating performance. [START_REF] Zhao | Understanding Electric Double-Layer Gating Based on Ionic Liquids: from Nanoscale to Macroscale[END_REF] In another work, the solvent effect on MoS 2 FET gated by a mixture of 1-butyl-3-methylimidazolium tetrafluoroborate and acetonitrile was studied, showing that the concentration of acetonitrile has a non-monotonic influence on the gating speed. [START_REF] Zhao | Adding Solvent into Ionic Liquid-Gated Transistor: The Anatomy of Enhanced Gating Performance[END_REF] The polarization effect is pronounced in ILs because of the strong local electric field. Although electrode polarization has been taken into account in recent simulation works on EDLs, the electrolytes are usually described by non-polarizable or even coarse-grained models, with only a few exceptions. [START_REF] Haskins | Evaluation of molecular dynamics simulation methods for ionic liquid electric double layers[END_REF][START_REF] Park | Interference of electrical double layers: Confinement effects on structure, dynamics, and screening of ionic liquids[END_REF] Recent advances in GPUaccelerated MD simulation make it possible to describe the electrified interface with more detailed models. [START_REF] Jiang | High-Performance Scalable Molecular Dynamics Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD[END_REF][START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF] In this work, full polarization is considered for both electrode and electrolytes by utilizing a constant-voltage simulation method together with the recently developed CL&Pol polarizable force field for ionic liquids. [START_REF] Goloviznina | A Transferable, Polarizable Force Field for Ionic Liquids[END_REF]

II. MODELS AND METHODOLOGIES

Constant-voltage simulations were performed to resemble the experimental conditions. That is, a voltage drop is pre-assigned between electrodes, while charge accumulated on the electrode is an observable, in contrast with the simulation condition where a constant charge density is pre-assigned to the electrodes and held fixed during the simulation. We note that since MoS 2 is a semiconductor, the constant charge density setup is likely to be more realistic than constant voltage method, which is strictly valid only for ideal conductor. However, recent works have shown that these two setups give similar interfacial structures and capacitances under low voltage drops, [START_REF] Haskins | Evaluation of molecular dynamics simulation methods for ionic liquid electric double layers[END_REF][START_REF] Park | Interference of electrical double layers: Confinement effects on structure, dynamics, and screening of ionic liquids[END_REF][START_REF] Wang | Evaluation of the constant potential method in simulating electric double-layer capacitors[END_REF] while a unique advantage of the constant voltage method is that it allows us to simulate the charging process more accurately. [START_REF] Merlet | Simulating Supercapacitors: Can We Model Electrodes As Constant Charge Surfaces?[END_REF] At the beginning of charging, a voltage drop is suddenly applied and charge starts accumulating at the electrode. The difficulty for constant-charge simulation during the charging process is that the correct charge density to be assigned on the electrode at each timestep is unknown.

The three-dimensional periodic image charge method recently proposed by Dwelle et al. [START_REF] Dwelle | Constant Potential, Electrochemically Active Boundary Conditions for Electrochemical Simulation[END_REF] was used in this work for constant-voltage simulations. Fig. 2 illustrates the geometry of the simulation box with ionic liquid confined between two MoS 2 plane electrodes. The electrodes are oriented along the xy plane and separated by d = 8.0 nm, which we measure as the distance between the S atoms in contact with ionic liquids on the cathode and the anode.

Considering the two electrodes as mirrors, each atomic charge in the electrolyte creates two oppositely signed image charges with respect to the two mirror planes. The image charges then generate higher-order image charges with respect to the mirrors, which can be captured by utilizing three-dimensional periodic boundary conditions with a length in the z direction of 2d = 16.0 nm. By applying an external electric field E ext on the electrolyte atoms, the voltage drop ∆V between two electrodes is guaranteed. 49

E ext = ∆V d (1)
Each electrode contains two layers of MoS 2 made of 13 × 8 unit cells. The length of simulation box in x and y directions are 4.109 and 4.380 nm, respectively. During image charge simulation, the electrode atoms are restrained at their original positions by a harmonic potential so that the locations of mirrors does not change.

U res = MoS2 i k x (x i -x i,0 ) 2 + k y (y i -y i,0 ) 2 + k z (z i -z i,0 ) 2
(2) herein, x i , y i and z i are the position of each Mo and S atom. x i,0 , y i,0 and z i,0 are the original position. k x , k y and k z are set to 0.1, 0.1 and 500.0 kcal mol -1 nm -2 , respectively. At each MD timestep, the positions of im-age charges are updated according to the coordinates of electrolyte atoms, immediately after the integration for positions.

The numbers of ion pairs for the systems studied are listed in Tab. I, which are determined based on the calculated densities of bulk ILs at 333 K and 1 bar. The simulation lengths are also listed in the table. Note that a trajectory of hundreds of nanoseconds is required to attain converged interfacial structures for these systems. The ILs are modeled with the recently developed Drude polarizable force field CL&Pol [START_REF] Goloviznina | A Transferable, Polarizable Force Field for Ionic Liquids[END_REF][START_REF] Pádua | Resolving dispersion and induction components for polarisable molecular simulations of ionic liquids[END_REF] with several refinements, whereas the electrode atoms are modeled as charge neutral particles with only LJ interactions. The details of force field and simulation conditions are given in section V. Near the electrodes, there are strong attractions between Drude particles and their images, which cause instability if a Drude particle comes too close to the electrode. Therefore, a wall potential is applied to all Drude particles of ILs to avoid this situation:

U wall =        Drude i 4 σ zi-z+ 12 - σ zi-z+ 6 + 0.25 , z i -z + < 1.1225σ Drude i 4 σ z--zi 12 - σ z--zi 6 + 0.25 , z --z i < 1.1225σ (3)
herein, and σ determine the strength and range of the wall potential, which are set to 0.5 kcal mol -1 and 0.15 nm, respectively. The z + and z -coordinates are the locations of cathode and anode, which are at 0 nm and 8.0 nm, respectively. We note that the wall potential is simply a safeguard. It decays to zero when the Drude particle is beyond 0.1684 nm from electrodes. The wall potential is only effective during the initial stage of the simulation when the interfacial strucure is evolving dramatically. By inspecting the converged trajectories, the minimum distance between Drude particles and electrode is 0.18 nm because of the vdW repulsion between heavy atoms. Therefore, this wall potential has no impact on the equilibrated trajectories.

III. RESULTS AND DISCUSSION

A. Structure of the interfaces Fig. 3 shows the distributions of centers of mass (COM) of imidazolium heads of [Im21] cations and [DCA] anions along z direction under voltage drops of 0 and 2 V. In the absence of voltage drop, [Im21][DCA] forms layered structures at the interface which oscillate for around 1.5 nm before reaching the bulk liquid. The symmetry of this distribution profile confirms the convergence of the trajectories. Herein, we use the peaks of the distribution of [DCA] to identify the location of layers at the interface, because the imidazolium head group can adopt different orientations near the electrodes. Three layers can be identified, which are located at 0.34, 0.72 and 1.09 nm from the electrodes. The [Im21] are more likely than the [DCA] to be found in the first layer. As the voltage increases to 2 V, the [Im21] and [DCA] ions then con- centrate near the anode and cathode, respectively. However, the locations of the three interfacial layers are not affected. The angles formed between ions and the plane of electrodes are analyzed to reveal the orientations of ions located in the first layer. Fig. 4 illustrates the definitions of three kinds of angles: the angle between the plane of an imidazolium ring and the electrode (α); the angle between the end-to-end vector of an octyl chain and the electrode (β); and the angle between the plane of a [DCA] anion and the electrode (γ). In [ImB1][DCA], the angle between the plane of benzyl ring and the electrode is also considered. The value of these angles range from 0°to 90°. The former corresponds to a parallel orientation to the electrode and the latter to a perpendicular one. is not homogeneous. There are polar-group rich domains and non-polar group rich domains, alternatively. This nano-scale segregation into polar and non-polar domains is well documented for longer-chain ionic liquids. [START_REF] Lopes | Nanostructural organization in ionic liquids[END_REF] Because the distance between two electrodes is in the order of nanometers, a confinement effect might have a contribution to this heterogeneity. [START_REF] Park | Interference of electrical double layers: Confinement effects on structure, dynamics, and screening of ionic liquids[END_REF] However, it should not have a pronounced effect on the structure of the EDLs. Because imidazolium rings and [DCA] always stay in the mutual vicinity in the bulk region, the polar and nonpolar segregation does not induce charge oscillations in bulk region. The applied voltage drop of 2 V promotes the enrichment of [Im81] and [DCA] near the anode and the cathode, respectively. The interfacial structure at the anode is similar to that at 0 V, except for the intensity of each peak. At the cathode, however, significant structural changes are observed. The interface with the cathode can not be simply identified as three layers any more. The octyl tails are pushed away from the first layer, therefore forming a non-polar rich region at around 0.8 nm from the cathode, which in turn induces strong oscillation of polar and non-polar domains.

The distributions of angles formed between the imidazolium head group, the octyl tail in at the cathode orient slightly more parallel to the electrode than those at the anode, which can be attributed to the electrostatic attraction between the three negatively charged nitrogen atoms in [DCA] and the cathode. However, the orientation of octyl tails is significantly altered by the applied voltage. It is shown that octyl tails at the cathode are more likely to be perpendicular to the interface, while the tails at the anode are still preferentially parallel to the electrode. This orientational change is counter-intuitive at first sight. Considering the attraction between the anode and the imidazolium head groups, it would be expected that the octyl tail leave the first layer near the anode so that more space is available for imidazolium groups. The orientational change observed can be attributed to the more favorable interactions between counter-ions than the interaction between alkyl tails and [DCA] anions. At the cathode, [DCA] occupy the first layer, then either imidazolium rings or octyl tails will be excluded from the first layer.

Here the more favorable interactions between oppositecharged ions results in the octyl chains being excluded.

A recent simulation work on butyltrimethylammonium bis(trifluoromethanesulfonyl)amide at a graphene elec-trode revealed a similar scenario, with butyl tails more likely to be found perpendicular to the electrode at the cathode. [START_REF] Katakura | Effect of Switching the Length of Alkyl Chains on Electric Double Layer Structure and Differential Capacitance at the Electrode Interface of Quaternary Ammonium-Based Ionic Liquids Studied Using Molecular Dynamics Simulation[END_REF] It was shown that the turning point in voltage for the butyl tails to be more likely perpendicular to the electrode is 1 V at the cathode, whereas the corresponding value is -5 V at anode. 7. Distributions of angles between imidazolium head groups, end-to-end vectors in octyl tails or [DCA] and xy plane in the [Im81][DCA] system under 0 V (upper plot) and 2 V (lower plot). The number inside the bracket denotes the z coordinate (in nm) of corresponding peak in density distribution, while the + anddenote the octyl tails bonded to the imidazolium head groups near the cathode and anode, respectively.

0
Fig. 8 shows the distributions of COM of imidazolium head groups and benzyl side groups in [ImB1] ions, and COM of [DCA] ions, along the z direction. In the absence of voltage drop, the first layer is mostly occupied by imidazolium rings and [DCA] anions. Only a small amount of benzyl groups are present in the first layer. The peak located at around 0.46 nm represents the benzyl groups that are attached to the imidazolium rings of the first layer. Also note that there is a peak for imidazolium head groups located at around 0.58 nm, which represents mainly the imidazolium rings attached to the benzyl groups in the first layer.

From the orientational distributions as shown in Fig. 9, both the imidazolium and benzyl groups located at 0.34 nm are parallel to the electrode surface, while the imidazolium rings located at 0.58 nm and the benzyl groups located at 0.46 nm are perpendicular to the electrodes. An applied voltage drop of 2 V has no effect to the shape of the orientational distributions, even thought the intensities of the density peaks are affected.

The [ImB1] cation adopts two main conformations in the first interfacial layer, as shown in Fig. 10. The first consists of imidazolium rings parallel to the electrode, while the benzyl ring is perpendicular with one hydrogen atom of the benzyl ring pointing toward the electrode. The second conformation corresponds to the benzyl ring parallel to the electrode, while the imidazolium ring is perpendicular to the electrode, with the methyl group pointing to the bulk electrolyte. In the absence of a voltage drop, the first conformation is dominant. Applying of voltage drop of 2 V induces exchange between these two conformations. In the first layer near the cathode, most of the [ImB1] that adopted the first conformation reorient into the second, because of repulsion between cathode and positively-charged imidazolium head-groups. At anode, the second conformation disappears.

B. Electrical properties of the EDLs

We are particularly interested in the link between the interfacial ordering and the electrical properties. Herein, the voltage profiles are extracted from the simulation trajectories and enable evaluation of capacitance. The volt- age profile along the z direction can be calculated by double integration of the charge distribution of the electrolytes, ρ(z), and of the charge density on the electrodes, Q, based on the one-dimensional Poisson equation:

V (z) = - 1 0 ρ(z)d 2 z + Qd ( 4
)
where 0 is the vacuum permittivity. Fig. 12 shows the charge distributions in the electrolytes. On the anode side, the first peak of positive charge is 0.27 nm away from the electrode. This distance is the same for all four electrolytes, because it is mainly determined by the vdW repulsion between S atoms of the electrode and the imidazolium ring. On the cathode side, the first peak of negative charge is 0.29 nm away from the electrode. However, we notice that a peak of positive charge is present at 0.25 nm away from the cathode, which is mainly due to hydrogen atoms from imidazolium cations. Even though the first layer at the cathode is dominated by [DCA] anions, hydrogen atoms from cations can approach the cathode more closely because of their smaller vdW radius, compared to those of heavy atoms.

Unlike in the fluctuating charge method, [START_REF] Siepmann | Influence of surface topology and electrostatic potential on water/electrode systems[END_REF] the charges on the electrodes are not directly obtainable in the present 3D periodic image-charge method, because the electrode atoms carry no charge by themselves. The effects of electrode polarization are described instead by the image charges and their infinite periodic images, together with the applied electric field. Herein, Q is evaluated from the positions of electrolyte atoms, 53

Q = n i=1 q i z i dA + ∆V 0 d (5)
where q i and z i are the charge and the z coordinate of each electrolyte atom, and A is the area of electrode.

Fig. 13 shows the voltage profiles along z across the electrolyte region with different ILs. The computed voltage drops between the electrodes equals to the preset values, and this validates the image charge method. The voltage profiles reach a plateau in the bulk electrolyte beyond ca. 2 nm away from the electrodes, showing the strong electrostatic screening of the EDLs. From the calculated charge densities on the electrodes and the voltage profiles, the integral capacitances of the EDLs at each electrode and the cell are evaluated as

C ± = Q V ± -V ± 0 (6) C cell = Q ∆V (7)
where V ± is the voltage drop between each electrode and the bulk region, and V ± 0 is the value of V ± in the absence of voltage drop between the electrodes. The calculated capacitances for different system are listed in Table II. For all ILs, the interfacial capacitance at the cathode is slightly larger than that at the anode. Among the four electrolytes, [ImB1][DCA] shows the largest cell capacitance, while [Im81][DCA] shows the smallest, the difference between them being quite small, only 8.5 %. The values of capacitance do not correlate with the molecular volumes of the ions, which indicates that the available space for counter-ions in the interfacial region is not the dominant factor for capacitance under voltage drop of 2 V. An applied voltage drop drives the cations and anions aggregating nearby anode and cathode, respectively. In order to inspect the charging mechanism, Fig. 14 shows the increase in the number of different species at the first layer near each electrode when voltage drop of 2 V is applied. An ion or fragment is considered to be located in the first layer if its COM is less than 0.5 nm away from the electrode.

Electrolyte C + (µF cm -2 ) C -(µF cm -2 ) C cell (µF cm -2 ) [Im21][DCA] 4.
In all of the systems, it is shown that the adsorption or desorption of [DCA] anions at the cathode or the anode dominates the behavior of ions at the interface, which can be attributed to the smaller size and therefore larger charge density of the [DCA] anion compared to imidazolium head-group. For [Im21], [ImB1] and [ImO1], the applied voltage drives imidazolium groups move into, or out of, the first layers of cathode or anode, which indicates a charing mechanism of exchange between counterions. However, the situation is different for [Im81]: at the cathode, the imidazolium head does not move out of the first layer. Instead, the tails are pushed out of the first layer, and more imidazolium head-groups enter into the first layer to occupy the created vacancies together with anions. At the anode, even though both imidazolium head groups and octyl chains become enriched in the first layer, the increment for octyl tails is larger than for imidazolium head-groups. Therefore, the charging mechanism for [Im81][DCA] is mainly the exchange between anions and octyl tails.

Another phenomenon worth noting is that the increase in the number of imidazolium head-groups at the cathode, and its decrease at the anode, is more significant for [ImB1] than for the other cations. This can be attributed to the orientational preference of [ImB1] cation, as discussed in previous section. As shown in Fig. 10, the conformation (b) dominates the first layer at the cathode. The benzyl group tends to be aligned parallel to the electrode and drives the imidazolium ring out of the first layer. At anode the conformation (a) dominates, with the benzyl group tending to be placed perpendicular to the electrode and occupying less space. The presence of the benzyl group promotes the separation of imidazolium head-groups and [DCA] anions at the first interfacial layer. As a result, [ImB1][DCA] is able to deliver slightly larger capacitance than the other ILs.

C. Role of the alkyl tail of [Im81]

It was recently reported that surface-active ILs promote enhanced capacitance by self-assembling into lamellar structures at the interface. [START_REF] Mao | Self-assembled nanostructures in ionic liquids facilitate charge storage at electrified interfaces[END_REF] In our simulation, even though polar and non-polar domain separation is observed in the bulk region, the [Im81] does not form pronounced layered structures near the anode. Instead, the octyl side chain is likely to occupy the first layer at the anode. Among the four ILs studied here, [Im81] delivers the smallest interfacial capacitance at the anode, as shown in Tab. II. It is therefore interesting to see whether the interfacial capacitance at the anode will increase if the octyl tails are forced to leave the first layer, so that more space is made available for imidazolium head-groups. We defined a reaction coordinate as the number of octyl tails in the first layer at the anode (N tail ). In order to make N tail differentiable with respect to atomic positions, the contribution of an atom i to its number density at position r is defined as a Gaussian function:

ρ i (r) = 1 σ √ 2π exp -|r -r i | 2σ 2 2 (8)
where r i is the position of atom i, and σ determines the radius of the atom, which we set to 0.05 nm. Then N tail is evaluated as

N tail = 1 2 C3-C8 i 1 6 erfc z 0 -z i √ 2σ (9)
where C3-C8 means the third to eighth carbon atoms in the octyl tail, z i is the z coordinate of each atom, and z 0 is the lower bound of the first layer at the anode, which we set to 7.5 nm based on the molecular distributions discussed in section III A. In order to force the octyl tail to leave the first layer, a harmonic biasing potential is applied to restrain the N tail to a target value:

U umbrella = 1 2 k(N tail -N target ) 2 (10)
where k is the harmonic force constant, which we set to 100.0 kJ mol -1 . Three simulations were performed with N target set to 0.0, 10.0 and 20.0, under the voltage drop of 2 V. Fig. 15 shows the number of imidazolium head-groups and [DCA] anions located in the first layer of the anode, when different N target are applied. The integral capacitance of the cell is also plotted. When N target is decreased, the octyl tails are forced to leave the first layer. As expected, more imidazolium head-groups are then able to occupy the first layer because of the free space created. However, it is found that the number of [DCA] ions in the first layer increases almost in parallel with the number of imidazolium groups. This indicates that the imidazolium groups and the [DCA] anions move into or out of the first layer in a correlated manner. Therefore, the charge that can be stored in the interface is not affected, which results in an essentially similar capacitance independently of the value of N target .

D. Charging performance

The speed of charge accumulation on the electrodes is one of the key properties of EDL devices, since it determines the power of EDLC and the gating speed of the FET. In order to inspect the charging process of the cell, initial configurations were taken from the system equilibrated under 0 V. Then a voltage drop of 2 V was applied to the electrodes and the evolution of charge density accumulated on the electrodes was extracted according to equation (5). performance for [ImO1] is close to that for [Im21] and is significantly better than those for [ImB1] and [Im81]. These results are surprising given that [ImO1] is three times more viscous than [Im21] and has comparable viscosity to [ImB1] and [Im81]. As such, we conclude that the charging performance is not strongly correlated with macroscopic dynamics. We interpret this results as being due to the charging processes of [Im21] involving only the displacement of small ions in adjacent layers, whereas that for [Im81] and [ImB1] are accompanied by conformational changes of cations at the interface. During the charging process, the octyl tails of [Im81] near the cathode are pushed away from the first layer and reorient to become perpendicular to the interface. For [ImB1], however, most of the cations of the first layer near the cathode change their conformations from (a) to (b), as shown in Fig. 10. These local conformational change slow down the charging of the EDL. In [ImO1], the 2-hydroxyethyl groups also change their orientations near anode, but the barrier for this conformational change is low because of the small size and flexibility. Therefore, the charging performance of [ImO1] is close to that of [Im21] and much higher than those of [Im81] and [ImB1], even though [ImO1] is much more viscous than [Im21]. From our simulations, the full charging of the cells happens at the time scale of nanoseconds. It is much shorter than experimental charing times reported for a IL-EDLC [START_REF] Chmiola | Anomalous increase in carbon capacitance at pore sizes less than 1 nanometer[END_REF] and a IL-FET, [START_REF] Fujimoto | Ionic-liquid component dependence of carrier injection and mobility for electric-double-layer organic thin-film transistors[END_REF] which are in the order of seconds. The discrepancy can be attributes to two factors: Firstly, in both the EDLC and FET mentioned above, porous electrodes are used, in which the diffusion of ions into or out of the electrodes dominates the charging speed; Secondly, the simulations were performed at 333 K in this work, which allows for much faster dynamics than at ambient temperature.

IV. CONCLUSIONS

We have performed polarizable MD simulations to study the interfaces between four imidazolium dicyanamide ionic liquid electrolytes and MoS 2 electrodes. The effect of side chain modifications to the structural, electrical and dynamic properties of the electric double layers are studied by using ethyl, octyl, benzyl and 2-hydroxyethyl side groups on the imidazolium cations ([Im21], [Im81], [ImB1] and [ImO1], respectively).

The octyl and benzyl groups are found to have prominent effects on the interfacial structure, while the 2hydroxyethyl group has almost no effect. In the absence of a voltage drop, the octyl tails of the [Im81] cation tend to lay parallel to the electrode in the first interfacial layer. An applied positive voltage drives the octyl tails out from the first layer, and they become perpendicular to the electrode, whereas an applied negative voltage has minor effect on the orientation of the octyl chains. The [ImB1] cation displays two major conformations in the first layer of the interface, one with the imidazolium head-group parallel to the electrode and the other with the benzyl group parallel to the electrode. In the absence of a voltage drop, the first conformation dominates. An applied positive voltage promotes a conformational change from the first to the second one.

The difference in interfacial structures results in distinct charging mechanisms. The charging mechanisms for [Im21], [ImB1] and [ImO1] consist mainly in exchange between counter-ions, whereas the mechanism for [Im81] is mainly the exchange between [DCA] anions and octyl tails. The separation between counter-ions is more pronounced for [ImB1][DCA] than for other ILs. This is due to the orientational preference of the [ImB1] cation. At cathode, the benzyl group lays parallel to the electrode, driving imidazolium group away from the first layer. While at anode, the benzyl group prefers to be perpendicular to the electrode, allowing imidazolium group to enter the first layer more easily. The benzyl group in [ImB1] assists with the separation of imidazolium group and [DCA] anion, therefore allows [ImB1][DCA] to deliver slightly larger capacitance.

In spite of the difference in molecular volumes and interfacial structures, all the four ILs deliver similar capacitances. It is found that the space available at the interface has little effect on the capacitance under a voltage drop of 2 V. Instead, the capacitance is mainly determined by the association between counter-ions. In [Im81], even when the octyl tails are forced out of the first layer, the space s created at the interface are occupied by neutral clusters of imidazolium head-groups and [DCA] anions, and as such the capacitance is not sharply affected. This result is comparable to the observation that the addition of organic solvents does not have significant effect on the interfacial capacitance, because of the balance between ionic concentration and the separability of the counter-ions. [START_REF] Burt | Capacitance of Nanoporous Carbon-Based Supercapacitors Is a Trade-Off between the Concentration and the Separability of the Ions[END_REF] From the analysis of charge accumulation on the electrodes during charging, it is shown that the charging speed with [ImO1] is close to that of [Im21], and is significantly larger than those of [Im81] and [ImB1], even thought the viscosity of [ImO1] is closer to that of [Im81] and [ImB1], than to [Im21], which is considerably lower. This difference in charging is attributed to the slow conformational changes of [Im81] and [ImB1] near the cathode during charging, and indicates that the charging performance does not correlate strongly with macroscopic dynamics. Instead, it is mainly determined by local reorientations and displacements at the interface.

This work emphasizes the importance of how microscopic molecular structures of ILs affect the properties of EDLs. In order to optimize the interfacial capacitance at a planar electrode, efforts should be undertaken on how to effectively suppress the association between counterions at the interface. In contrast, the capacitance is largely insensitive to the molecular volume for ILs of the same family. The present results also point out that bulk dynamic properties like viscosity do not determine by themselves the charging performance of EDLs. Instead, the microscopic conformational changes happening at the interface assume an important role in the charging process, and should be taken into account in the rational design of IL electrolytes. In the present, all cations share the same imidazolium head. It is therefore intriguing to see whether these conclusions still stand for ILs belonging to different families (e.g. tetraalkylammonium, phosphonium), considering that the association strength between counter-ions can be very different.

V. FORCE FIELD AND SIMULATION DETAILS

The van der Waals (vdW) parameters for MoS 2 were obtained from the work of Sresht et. al. [START_REF] Sresht | Quantitative Modeling of MoS2-Solvent Interfaces: Predicting Contact Angles and Exfoliation Performance using Molecular Dynamics[END_REF] Since the image-charge method is used to describe the Coulomb interactions between MoS 2 and the electrolytes, the electrode atoms carry no charge. The recently-developed CL&Pol 46,50 polarizable force field is used to describe the ILs. In order to avoid the polarization catastrophe, the electrostatic interactions between hydrogen atoms in hydroxyl group and all Drude dipoles are damped by using a Tang-Tonnies function. [START_REF] Goloviznina | Transferable, polarizable force field for electrolytes, protic ionic liquids and deep eutectic solvents[END_REF] The Thole screening parameter and non-bonded parameters for hydroxyl and benzyl group are refined to provide better agreement with experimental viscosity data.

The CL&Pol model was developed on top of the nonpolarizable CL&P force field 58 by introducing Drude induced dipoles to heavy atoms. After the introduction of Drude induced dipoles, the Lennard-Jones (LJ) parameters are scaled, because the non-polarizable force field had been optimized to account for part of the polarization implicitly. A predictive scheme to evaluate these scaling factors for LJ well-depths from the charge and dipole moment of fragments was proposed during the development of CL&Pol. [START_REF] Goloviznina | A Transferable, Polarizable Force Field for Ionic Liquids[END_REF] The scaling factor for the LJ radius is empirically determined to be 0.985 based on the calculated densities.

The refined force field for ILs was validated by calculating the viscosity and density of of several ILs at 333 K and 1 bar. The viscosities are calculated with periodic perturbation method, [START_REF] Hess | Determining the shear viscosity of model liquids from molecular dynamics simulations[END_REF] which has proven to be accurate and efficient for viscous liquids like ionic liquids, or liquids under extreme pressure. [START_REF] Gong | Extension of TEAM Force-Field Database to Ionic Liquids[END_REF][START_REF] Gong | Pressure-viscosity relation of 2,2,4trimethylhexane predicted using all-atom TEAM force field[END_REF] In this method, an acceleration along the x direction, with periodicity in z, is applied in an otherwise normal N pT simulation:

a x (z) = A cos 2πz l z (11)
where l z is the z-length of the simulation box. In steady state, the acceleration will generate a velocity gradient:

v x (z) = W cos 2πz l z (12)
from which the velocity amplitude W can be calculated by ensemble averaging of the simulation data. Via Stokes' equation, the amplitude is related to the viscosity η by:

W = Aρ η l z 2π 2 (13)
Viscosities were calculated at four accelerations (0.01, 0.02, 0.03 and 0.04 mm ps -2 ), from which a linear extrapolation towards zero acceleration is performed to obtain the zero-shear viscosity. [START_REF] Gong | Extension of TEAM Force-Field Database to Ionic Liquids[END_REF][START_REF] Zhao | On the accuracy of predicting shear viscosity of molecular liquids using the periodic perturbation method[END_REF] The density and viscosity are calculated with different values of the Thole parameter to determine its effect. As shown in Fig. 17, the Thole parameter has little effect on these properties, with the exception that the density calculated with extended-Langrangian (EL) approach is smaller than the result of self-consistent field (SCF) method when the Thole parameter is set to the usual value of 2.6. By inspecting the root-mean-squared force (RMSF) on Drude particles after 20 SCF iterations at the last step, it is shown that decreasing the Thole parameter makes the SCF iterations more robust. Therefore, the Thole parameter is chosen to be 1.8 in this work to enhance the stability of the Drude model without affecting the thermodynamic and transport properties. In the CL&Pol force field, the LJ well-depth between cation-anion, cation-neutral, anion-neutral and neutralneutral fragments are scaled down using the predictive scheme, [START_REF] Goloviznina | A Transferable, Polarizable Force Field for Ionic Liquids[END_REF] whereas the LJ well-depth between fragments with charges of the same sign were kept untouched (since these are seldom in the respective vicinity). In this work, the LJ well-depth between fragments with charges of the same sign are also scaled down based on the predictive scheme. [START_REF] Goloviznina | A Transferable, Polarizable Force Field for Ionic Liquids[END_REF] However, the predictive scheme requires a minimal-energy distance between fragments optimized by QM calculation, which is not applicable for same-charge fragments. As such, we estimate the distance between same-charge fragments by d(frag1, frag2) = d(frag1, butane) + d(frag2, butane) -d(butane, butane) (14) After scaling down the LJ well-depth of same-charge fragments, the experimental viscosities for alkylimidazolium dicyanamide were more accurately reproduced, as shown in Table III.

The Drude polarizable trajectories were generated with an integrator using the extended-Langrangian approach, where the COM motion of real atoms and attached Drude particles are thermalized at 333 K, unless otherwise specified, while the relative motion between atoms and their Drude particles are thermalized at 1 K to approximate the SCF relaxation of Drude dipoles. [START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF] The electrode atoms are thermalized using a Langevin thermostat with a friction coefficient of 5 ps -1 . The IL atoms are thermalized using the recently developed temperature-grouped Nosé-Hoover (TGNH) thermostat in order to attain more reliable dynamic properties. [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models : Temperature Grouped Dual-Nose-Hoover Thermostat[END_REF] With the conventional dual-Nosé-Hoover (NH) thermostat, [START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF][START_REF] Martyna | Nosé-Hoover chains: The canonical ensemble via continuous dynamics[END_REF] the temperatures of different degrees of freedoms (DOFs) cannot be properly maintained. Because of the continuous heat flow from DOFs of real atoms to Drude particles, the temperature of molecular COM motions becomes significantly higher than that of other DOFs, which causes underestimation of density and overestimation of diffusion. As a remedy, the TGNH thermostat thermalizes the DOFs of molecular COM motions, atomic internal motions relative to molecular COM and Drude relative motions separately, which gives more accurate estimations of density and diffusion coefficient. [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models : Temperature Grouped Dual-Nose-Hoover Thermostat[END_REF] Three Nosé-Hoover chains were used for the TGNH thermostat. The coupling strength for real-atom DOFs and Drude DOFs are 10 ps -1 and 40 ps -1 , respectively. As shown in Table III, the conventional NH thermostat is not able to correctly thermalize Drude-polarizable systems. The temperature of molecular COM motion is over 50 K higher than the preset value, which results in significant underestimation of viscosity and density. In contrast, the TGNH thermostat is able to maintain the average molecular translational temperature, thus predicting density to be the same as the one obtained with the Langevin thermostat.

The vdW and the real space electrostatic interactions are truncated at 1.2 nm. The particle mesh Ewald method [START_REF] Darden | Particle mesh Ewald: An N•log(N) method for Ewald sums in large systems[END_REF] is applied to handle the long-range electrostatic interactions in reciprocal space. The integration time step is 1 fs. All bonds involving hydrogen and the C-O-H angles in [ImO1] are constrained by using SHAKE algorithm. [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF] In N pT simulations for determining density, viscosity and permittivity, a Monte Carlo barostat [START_REF] Chow | Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling[END_REF][START_REF] Åqvist | Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm[END_REF] was used to maintain pressure, wherein the volume of simulation box and positions of molecular centers were scaled every 100 ps. The uncertainties reported in this work were estimated as the standard errors of five blocks using the block-average method. The extended-Langrangian simulations were performed with OpenMM 7.4.1, 45,74 while the SCF simulations were performed using GROMACS 2016. [START_REF] Hess | GRO-MACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation[END_REF] OpenMM is extensively accelerated on GPU processors, allowing sufficiently long simulation trajectories to be generated. The methods used in this work were implemented as an OpenMM plug-in. The source code is hosted at https://github.com/z-gong/ openmm-velocityVerlet. The computational tools and data that support the findings are made public in the above-mentioned repositories and are available from the corresponding author upon reasonable request.

Zheng

  Gong and Agilio A. H. Padua a) Laboratoire de Chimie, École Normale Supérieure de Lyon & CNRS, 69364 Lyon, France (Dated: 28 January 2021)

  ][DCA]) electrolytes at planar MoS 2 electrodes using polarizable MD simulations. Four cations are considered: 1-ethyl-3-methylimidazolium ([Im21]), 1-octyl-3-methylimidazolium ([Im81]), 1-benzyl-3methylimidazolium ([ImB1]) and 1-(2-hydroxyethyl)-3methylimidazolium ([ImO1]), as depicted in Fig. 1. The ethyl, octyl, benzyl and 2-hydroxyethyl groups are representative for different shape, size and hydrophobicity.

FIG. 1 .

 1 FIG. 1. Structures of ions considered in this work.

FIG. 2 .

 2 FIG. 2. System geometry and setup for image charge simulations of ILs confined between two MoS 2 electrodes. The purple points are image charges. The blue box shows the periodic cell boundary. Each layer of MoS 2 is covalently bonded to its images in the x and y directions.

FIG. 3 .

 3 FIG. 3. Distributions of COM of different fragments along z direction in the [Im21][DCA] system under 0 V (upper plot) and 2 V (lower plot).

FIG. 4 .FIG. 5 .

 45 FIG. 4. Illustration of the angles formed between imidazolium head, octyl tail or [DCA] anion and the electrode. The vertical line at the left side represents the electrode in xy plane. Hydrogen atoms are omitted for clarity.

10 FIG. 6 .

 106 FIG.6. Distributions of COM of different fragments along the z direction in [Im81][DCA] system under 0 V (upper plot) and 2 V (lower plot).

10 FIG. 8 .

 108 FIG. 8. Distributions of COM of different fragments along the z direction in the [ImB1][DCA] system under 0 V (upper plot) and 2 V (lower plot).

Fig. 11 shows

 11 Fig. 11 shows the distributions of COM of imidazolium rings and oxygen atoms in [ImO1] ions, and also COM of [DCA] ions, along z. The distributions of imidazolium rings and [DCA] anions are almost identical to those in [Im21][DCA]. The oxygen atoms of the hydroxyl groups are more likely to stay in the first layer near the cathode.However, because of its flexibility, the location of the hydroxyl group has a minor effect on the structural arrangement of imidazolium head groups and anions. The orientations of ions are also similar to those in [Im21][DCA], and thus are not shown here.

FIG. 9 .

 9 FIG. 9. Distributions of angles between imidazolium headgroups or benzyl rings or [DCA], and the xy plane in the [ImB1][DCA] system under 0 V (upper plot) and 2 V (lower plot). The number inside the bracket denotes the z coordinate (in nm) of corresponding peak in density distribution. The profiles for ions near the anode are not shown for clarity.

FIG. 11 .

 11 FIG. 11. Distributions of COM of different fragments along the z direction in [ImO1][DCA] system under 0 V (upper plot) and 2 V (lower plot).

FIG. 12 .

 12 FIG.12. Charge density distributions along the z direction for different electrolytes under voltage drop of 2 V.

FIG. 13 .

 13 FIG.13. Voltage profiles along the z direction for different electrolytes under voltage drops of 0 V (dashed lines) and 2 V (solid lines).

FIG. 14 .

 14 FIG. 14.The increase in the number of different species in the first layer near cathode (filled bars) or anode (unfilled bars) when voltage drop of 2 V is applied. The tail of [Im21] is not considered.

Fig. 16 FIG. 15 .

 1615 FIG. 15. Number of imidazolium head-groups and [DCA] anions in the first layer near the anode, and integral capacitance of the cell at different settings of Ntarget for [Im81][DCA], under a voltage drop of 2 V.

FIG. 16 .

 16 FIG.16. Time evolution of the charge densities accumulated at the electrodes for different electrolytes when voltage drop of 2 V is applied. The horizontal dashed lines denote the converged values of charge densities.

  FIG. 17. Calculated density, viscosity, and the RMSF on Drude particles for bulk [Im21][DCA] at 333 K with different values of the Thole parameter. The uncertainty of density is too small to be shown.

TABLE I .

 I The number of ion pairs and trajectory length for each electrolyte.

	Electrolyte Number Length (ns)
	[Im21][DCA]	511	180
	[Im81][DCA]	312	480
	[ImB1][DCA]	392	360
	[ImO1][DCA] 503	240

TABLE II .

 II Calculated integral capacitances of EDLs under a voltage drop of 2 V.

TABLE III .

 III Calculated viscosity, density and molecular COM temperature for several ILs with the refined CL&Pol force field using different thermostats. The uncertainties for calculated density are smaller than 0.001 g cm -3 . The experimental data[START_REF] Neves | Systematic study of the thermophysical properties of imidazolium-based ionic liquids with cyano-functionalized anions[END_REF][START_REF] Karpińska | Separation of ethylbenzene/styrene systems using ionic liquids in ternary lle[END_REF][START_REF] Xu | Solubility of co2 in aqueous mixtures of monoethanolamine and dicyanamide-based ionic liquids[END_REF][START_REF] Paduszyński | An effect of cation functionalization on thermophysical properties of ionic liquids and solubility of glucose in them -measurements and pc-saft calculations[END_REF] are listed in the brackets for comparison.Ionic liquid Thermostat T (K) TCOM (K) Viscosity (mPa s) Density (g cm -3 )

	[Im21][DCA] TGNH	333	333.3±0.4 6.3±0.3 (6.95)	1.088 (1.081)
	Langevin	333	334.9±0.9 -	1.088 (1.081)
	NH	333	388.8±1.3 2.6±0.1 (6.95)	1.046 (1.081)
	[Im61][DCA] TGNH	333	334.5±0.3 17.0±1.1 (15.0) 0.996 (1.008)
	[Im81][DCA] TGNH	333	332.7±0.2 25.5±1.9	0.976
	[ImB1][DCA] TGNH	333	333.7±0.5 26.5±2.3	1.126
	Langevin	298	297.9±0.9 -	1.157 (1.158)
	[ImO1][DCA] TGNH	333	333.4±0.3 19.6±1.2 (20.7) 1.166 (1.164)
	DATA SHARING

ACKNOWLEDGMENTS

This work was supported through funding from the Agence Nationale de la Recherche (grant ANR-18-CE09-0018). A CC BY or equivalent licence is applied to Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions. Computational resources were provided by IDRIS (national computing center of the CNRS) and the PSMN (computing center at ENS de Lyon).


Laurent Mathieu
email: laurent_tom2@yahoo.fr

Soryapong Plang

Nicolas De L'escalopier

James Charles Murison

Christophe Gaillard

Antoine Bertani

Frédéric Rongieras

Soft tissue coverage using pedicled flap in combat zone: a case series

Keywords: Pedicled flaps, Limited resources, War surgery, Reconstruction, Training

Background: Soft tissue reconstruction is typically conducted after evacuation from theater of operations. If circumstances do not allow timely evacuation, however, defect site may need to be reconstructed in the combat zone. Case presentation: A total of 41 patients with extremity soft tissue defect were treated using pedicled flaps by a single orthopedic surgeon during four deployments in Chad, Afghanistan and Mali between 2010 and 2017. The mean age was 25.6 years. A total of 46 injury sites in extremities required flap coverage: 19 combat-related injuries (CRIs) and 27 non-combat related injuries (NCRIs). Twenty of the injury sites were infected. Overall, 63 pedicled flap transfers were carried out: 15 muscle flaps, 35 local fasciocutaneous flaps and 13 distant fasciocutaneous flaps. The flap types used did not differ for CRIs or NCRIs. Mean follow-up was 71 days. Complications included deep infection (n = 6), flap failure (n = 1) and partial flap necrosis (n = 1). Limb salvage rate was 92.7% (38/41).

Conclusions: Soft tissue defect can be managed with simple pedicled flaps in theatre of operations if needed. Basic reconstructive procedures should be part of the training for military orthopedic surgeons.

Background

Modern conflicts present military surgeons with a high volume of extremity injuries requiring flap coverage [START_REF] Sabino | A decade of conflict: flap coverage options in traumatic war-related extremity reconstruction[END_REF]. Various studies have evaluated soft tissue coverage outcomes in U.S. military personnel after evacuation from combat zones [START_REF] Sabino | A decade of conflict: flap coverage options in traumatic war-related extremity reconstruction[END_REF][START_REF] Connolly | Changing paradigms in lower extremity reconstruction in war-related injuries[END_REF]. These reconstructive procedures performed in specialized centers by multidisciplinary teams are rarely accessible to local, national conflicts since they cannot be evacuated from the theater of combat [START_REF] Klem | Microvascular reconstructive surgery in operations Iraqi and enduring freedom: the US military experience performing free flaps in a combat zone[END_REF].

Due to the lack of plastic surgeons in the battlefield, orthopedic surgeons often find themselves alone when managing local patients who present with complex extremity injuries requiring multi-tissue reconstruction. In this austere environment, they routinely perform soft tissue coverage procedures using "simple, reliable and replicable" techniques [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF][START_REF] Barbier | Could the orthopaedic surgeon deployed in auster setting perform flaps on the leg?[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF].

In this case series, we report the use of pedicled flap transfers in combat zone medical treatment facilities (MTFs).

Case presentation

Data collection and analysis

A total of 41 patients with extremity soft tissue defects were treated using pedicled flaps by a single orthopedic surgeon (LM) between 2010 and 2017 during four deployments as a member of forward surgical teams (FST) in Chad and Mali, and in a combat support hospital (CSH) in Afghanistan. Three types of pedicled flaps were used: muscle flaps, local fasciocutaneous flaps, and distant (fascio) cutaneous flaps.

Postoperative complications are presented using the Clavien-Dindo classification [START_REF] Clavien | The Clavien-Dindo classification of surgical complications: five-year experience[END_REF]. Outcome measures included flap loss, partial flap necrosis and early infections. Flap loss was defined as a need for coverage-revision surgery. Partial flap necrosis was defined as necrosis that necessitated surgical debridement but did not require additional coverage surgery. Early infections were defined as a wound infection at the coverage site within 2 weeks of flap transfer that required a return to the operating theater. Time to additional skin grafting on muscle flaps or the donor sites of fasciocutaneous flaps was also analyzed.

Endpoint assessment included limb/finger salvage and late complications due to infection. Injuries were categorized based on energy trauma and wound contamination to combat-related injuries (CRIs) versus non-combat related injuries (NCRIs). Categorical variables were compared using Fisher's exact test. P < 0.05 was considered statistically significant.

Baseline characteristics

The study included 41 subjects (35 male and 6 female patients, with a mean age of 25.6 ± 15.0 years). The origin of the injuries was non-ballistic trauma in 18 cases, ballistic trauma in 15 cases, osteomyelitis in 6 cases and burn in 2 cases. CRIs occurred significantly more frequently in Afghanistan than Chad and Mali (13/15 versus 2/15, P = 0.001). Four patients presented with multiple lesions. Thus, a total of 46 injury sites (19 CRIs and 27 NCRIs) required flap reconstruction. Most injuries were located in the legs (n = 23) or hands (n = 15), and with open fractures ((n = 27, Table 1). Twenty of the 46 injury sites were infected.

Pattern analysis of flap surgery

The average number of debridement procedures per injury was 1.8 ± 0.9. Serial debridement was required prior to flap coverage in 29 injuries; negative wound pressure therapy was used between debridement sessions in 12 injuries. Primary debridement and flap coverage were carried out simultaneously in 17 injuries.

Two locoregional flaps were combined in 10 large defects, and 9 simultaneous distant abdominal flaps were required to cover a burn injury of both hands in a single patient. Thus, a total of 63 pedicled flap transfers were performed: 15 muscle flaps, 35 local fasciocutaneous flaps and 13 distant fasciocutaneous flaps. The flap types chosen for treatment did not differ between CRIs and NCRIs. Muscle flaps were mostly used for proximal and mid-tibia coverage; distal tibia coverage was achieved by transposition and island fasciocutaneous flaps (Fig. 1). Three thumb reconstructions were carried out using island digital flaps. Distant flaps were exclusively used for hand and forearm coverage (Fig. 2). Additional skin grafting was required in 40 of the 63 flap transfers and performed after a mean delay of 5 ± 5 days. Associated procedures included external fixation (n = 21), internal fixation (n = 11), bone reconstruction using the induced membrane technique (n = 10) and tendon repair (n = 3). Secondary division of the distant flap pedicle was carried out after three weeks.

Treatment outcome

The mean follow-up time was 71 ± 95 days. There were nine type III Clavien-Dindo complications: one abdominal flap loss, one partial necrosis of a groin flap (Fig. 2a), one knee joint-fluid fistula (Fig. 1a) and six early infections. All other complications were type I: two minimal marginal losses of fasciocutaneous flaps and one dehiscence at the donor site of a groin flap.

Three patients in Mali were offered a late amputation due to severe persistent infection; one of them declined below knee amputation for septic non-union of the tibia as he was able to walk on crutches without pain. Three patients had a chronic pus-fistula under their muscle flap reconstruction related to tibia infection: two had been treated for a neglected open fracture and one for severe osteomyelitis. Limb salvage rate was 92.7% (38/ 41).

Discussion

Battlefield MTFs provide life-and limb-saving care through damage-control surgical procedures. The conditions in the field are typically ill-suited to support To our knowledge, this is the first published case series of soft tissue coverage by a single surgeon over multiple tours of duty in various locations. In our opinion, this feature reduces the impact of confounding factors on the results (inter-surgeon variability) and enhances the applicability of the results across multiple types of theater of operations.

In modern reconstructive units within wellprovisioned healthcare facilities, soft tissue reconstruction is individually tailored to the wound, available and reliable flap sources, associated injuries and specific rehabilitation goals for each patient [START_REF] Sabino | A decade of conflict: flap coverage options in traumatic war-related extremity reconstruction[END_REF]. In battlefield MTF environment, however, the choices for soft tissue coverage methods are restricted. Factors that affect treatment decisions include the surgeon's expertise and available resources (e.g., surgical equipment, antibiotics, laboratory analysis facilities and the number of available beds) [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF]. For these reasons, the simplest solution for coverage is always preferred [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF]. In our experience, pedicled flaps combined with skin grafts allowed reconstruction of almost all types soft tissue extremity injuries, even large ones. Additional skin grafts were carried out together with flap transfers in most patients but were deferred in cases at risk of early infection. Simultaneous local, or distant, pedicled flaps were used successfully in this cohort as an alternative to free transfers [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF]. This decision was based on our limited experience with such procedures, as well as the following 2 considerations. First, free flaps require microsurgical techniques and specific post-operative care, and thus could hardly be carried out by non-specialized healthcare teams. Second, free flap surgery requires an extended period in theater that can jeopardize the operational activity of a forward surgical facility [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF].

In the present study, the failure rate was low for both CRIs and NCRIs. The overall success rate for flap coverage performed in the field was over 90%, a figure comparable to various other authors' reports about warrelated extremity reconstructions performed in patients after evacuation out of the combat zone [START_REF] Sabino | A decade of conflict: flap coverage options in traumatic war-related extremity reconstruction[END_REF][START_REF] Connolly | Changing paradigms in lower extremity reconstruction in war-related injuries[END_REF]. In our opinion, the relatively high success rate could be attributed to the near-exclusive use of simple, reliable, and replicable pedicled flap transfers, which were perfectly suited to both CRIs and NCRIs in a non-specialist surgical MTF [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF][START_REF] Barbier | Could the orthopaedic surgeon deployed in auster setting perform flaps on the leg?[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF]. Transposition, or rotational, fasciocutaneous flaps and muscle flaps were the two types most often used, regardless of the injury location or cause. Since these flaps do not require pedicle dissection, the procedure could be readily performed by surgeons with no specialized training in plastic surgery.

Other flap types (e.g., perforator and distant flaps) were also used in this study. Perforator flaps are not recommended for CRI treatment due to the extensive soft tissue injury and the potential violation of fascial planes and perforators caused by the projectile's kinetic energy [START_REF] Connolly | Changing paradigms in lower extremity reconstruction in war-related injuries[END_REF]. They were mostly employed to treat NCRIs; they were dissected as island flaps with a large adipofascial pedicle (designed with a ratio L/l < 4) following Doppler examination. When using such flaps, non-specialized surgeon should be aware that the flap design is based on the location of the vascular territory of the perforator as well as the perforator flow direction [START_REF] Saint-Cyr | The perforasome theory: vascular anatomy and clinical implications[END_REF]. Propeller perforator flaps were never used because they were too technically demanding [START_REF] Sisti | Propeller flaps: a literature review[END_REF]. Distant flaps were indicated to salvage upper limbs, mostly at the hand level, even though they required a minimum 3-week hospital stay. By contrast, we did not use cross-leg flaps for tibia coverage as our experience has led us to conclude that an amputation should be considered in absence of available local flaps in such austere healthcare settings [START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF].

There were 12 flap complications, half of which were due to early infections, highlighting the importance of infection control in both CRIs and NCRIs. Regardless of the method of soft tissue reconstruction, adequate debridement of necrotic or infected tissue is critical for the overall success of any reconstructive modality [START_REF] Connolly | Changing paradigms in lower extremity reconstruction in war-related injuries[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF]. Highly contaminated war wounds usually require serial debridement and broad-spectrum intravenous antibiotics prior to definitive coverage [START_REF] Connolly | Changing paradigms in lower extremity reconstruction in war-related injuries[END_REF]. Management of NCRIs, such as neglected open fractures, septic non-unions, and chronic osteomyelitis, also follow the same therapeutic rules. Within the limitations of MTF conditions, therefore, treatment of bone infections should be undertaken with caution if extended courses of antibiotics, sequential procedures, and close monitoring (for several months) will not be possible [START_REF] Mathieu | Soft-tissue coverage of war extremity injuries: the use of pedicle flap transfers in a combat support hospital[END_REF][START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF].

The current study demonstrated that an orthopedic surgeon with basic knowledge in local vascular anatomy is able to harvest an appropriate local, regional or distant pedicle flap to manage the majority of soft tissue defects in the field. A supplemental training in reconstructive techniques would, however, be required for deployed orthopedic surgeons who are not familiar with extremity reconstruction. In the French Military Health Service such training is now included in the Advanced Course for Deployment Surgery [START_REF] Bonnet | The French Advanced Course for Deployment Surgery (ACDS) called Cours Avancé de Chirurgie en Mission Extérieure (CACHIRMEX): history of its development and future prospects[END_REF]. During the training, basic pedicled flap transfers and bone reconstruction techniques that can be performed in healthcare facilities with limited resources are learned through lectures, hands-on exercises on cadavers and case studies [START_REF] Mathieu | Type III open tibia fractures in low resources setting. Part 2: soft-tissue coverage with simple, reliable and replicable pedicle flaps[END_REF].

This study has several limitations. First, the study population was heterogeneous in terms of factors like injury mechanism and the time elapsed before management. Second, indications for the different pedicle flaps are open for discussion as they reflected only the views and experience of one surgeon. Third, the short followup made it impossible to evaluate long-term limb salvage, bone infection control and achievement of bone union.

Conclusions

Pedicled flap transfers are safe and useful procedures suitable for soft tissue coverage within forward surgical units. All military orthopedic surgeons should be trained to perform such basic reconstruction techniques. Except perhaps in cases of pre-existing bone infection, these techniques permit limb salvage in most open extremity soft tissue injuries encountered in the field.

3 ) 0

 30 Ankle and foot (n = 3) 1 (33.3) 0 2 (66.7) 0 Elbow and forearm (n = 5) 3 (60.0) 1 (20.0) 1 (20.0) 0 Hand (n = 15) 5 (33.3) 1 (6.7) 6 (40.0) 3 (20.0) Total (n = 46) 27 (58.7) 6 (13.0) 10 (21.7) 3 (6.5) specialized plastic surgery. Therefore, studies reporting the results of flap reconstruction in theatres of operations are rare. Barbier et al. [5] analyzed the use of muscle and rotational fasciocutaneous flaps for recent, or neglected, open-tibia fractures, septic nonunions, and osteomyelitis. We have also previously reported a series of 35 pedicled flap transfers for soft tissue reconstruction of various CRIs in Afghanistan [4]. Klem et al. [3] reported on a case series of microsurgical free-flap procedures, performed by plastic or ear, nose and throat surgeons together with orthopedic surgeons, in U.S. CSHs that were deployed in Iraq and Afghanistan.

Fig. 1

 1 Fig. 1 Locoregional flap transfers in the lower extremity. a Simultaneous lateral and medial gastrocnemius flaps for a blast injury of the knee. b Distally based great saphenous flap to cover a distal tibia open fracture in a 10-year-old patient

Table 1

 1 Injury pattern and distribution [n(%)]

	Injury sites	Open fracture	Osteomyelitis	Soft tissue injury	Burn injury
	Knee and leg (n = 23)	18 (78.3)	4 (17.4)	1 (4.

Acknowledgments

The authors are grateful to Jennifer Dandrea Palethorpe for assistance in preparing this manuscript, and to all medical and paramedical staff of the various MTFs who took part in the treatment of the patients.

Funding

Not applicable.

Availability of data and materials

The datasets analyzed during the current study are available from the corresponding author upon reasonable request.

Ethics approval and consent to participate

All patients gave oral consent to participate in this case study. This study was approved by the Department of Surgery of the French Military Health Service Academy (2019-0901-001).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.