The role of solvent interfacial structural ordering in maintaining stable graphene dispersions

Liquid phase exfoliation is the most promising method for the low-cost, scalable production of two-dimensional nanosheets from their bulk counterparts. Extensive exfoliation occurs in most solvents due to the huge amount of energy introduced by sonication or shear mixing. However, the subsequent dispersion is not always stable, with extensive reaggregation occurring in some solvents. Identifying the optimal solvent for a particular layered material is difficult and requires a fundamental understanding of the mechanism involved in maintaining a stable dispersion. Here, we use molecular dynamics calculations to show that when graphene is immersed in a solvent, distinct solvation layers are formed irrespective of the choice of solvent and their formation is energetically favourable for all considered solvents. However, energetic considerations such as these do not explain the experimental solvent-dependence of the dispersion concentration. Instead, we find that solvents with high diffusion coefficients parallel to the graphene layer result in the lowest experimental concentration of graphene in solution. This can be explained by the enhanced ease of reaggregation in these solvents. Solvents with smaller diffusion coefficients result in higher experimental graphene concentrations as reaggregation is prevented. In the low diffusion limit, however, this relationship breaks down. We suggest that here the concentration of graphene in solution depends primarily on the separation efficiency of the initial exfoliation step. Based on this, we predict that the concentration of exfoliated graphene in solvents such as benzaldehyde and quinoline, which have low diffusion constants, can be increased dramatically by careful tuning of the experimental sonication parameters.


Introduction
The liquid-phase exfoliation (LPE) of layered materials is a scalable top-down method to produce industrial quantities of monolayer or few-layer two-dimensional (2D) nanosheets at a reasonable cost [1][2][3][4].Solventdispersed 2D nanosheets can find applications as inks for printed flexible electronics [5,6], energy generation and storage [7], sensing [8] and fire resistant coatings [9].The ability of LPE to produce large quantities compensates for relatively small lateral size flakes and low yield achieved compared to other exfoliation methods [10].However, the quality and quantity of the exfoliated layers depends dramatically on the solvent used [11,12].A careful choice of both the solvent and the sonication parameters can improve the monolayer yield and increase the lateral dimensions achieved, but a fundamental lack of understanding of the mechanism driving the exfoliation and stabilization processes impedes this optimisation procedure [13,14].
These two processes, exfoliation and stabilization, can be considered separately even though they occur effectively simultaneously experimentally.Exfoliation occurs when the van der Waals interactions binding the layers together are overcome via large shear forces that are introduced in the presence of a solvent through either sonication, high-shear mixing or wetball milling [10].The solvent has been suggested to play a relatively minor role in this initial separation of layers, with extensive exfoliation occurring in most solvents due to the huge amounts of energy involved [15,16].The solvent then stabilizes the exfoliated nanosheets preventing reaggregation and sedimentation.
The choice of solvent now plays a critical role.Solvent exchange methods, where different solvents are used in the exfoliation and stabilization steps, find that those solvents which do not perform well in the initial exfoliation step may perform very well in maintaining a stable dispersion of already exfoliated monolayers [17][18][19].Evidently, the solvent plays a different role in the two steps but as they occur effectively simultaneously, their exact role is difficult to disentangle.
The majority of investigations in this regard concern the exfoliation of graphite.Hernandez et al. compared the performance of a wide selection of solvents and found that the maximum graphene concentration was achieved by exfoliation in cyclopentanone, resulting in 8.5 ± 1.2 µg/ml [20].Exfoliation in non-polar solvents, including toluene, heptane, hexane and pentane, resulted in much lower graphene concentrations of 0.8 ± 0.4 µg/ml, 0.3 ± 0.4 µg/ml, 0.2 ± 0.1 µg/ml and 0.16 ± 0.05 µg/ml, respectively.However, some highly polar molecules such as water and formamide also perform poorly, with the best performing solvents having a slightly polar nature (although not all slightly polar solvents perform well).The effectiveness of LPE for other layered materials, such as WS 2 , MoS 2 and h-BN, is also strongly dependent on the choice of solvent [2,21].
Clearly, the inability to identify the optimal solvent for the exfoliation of a particular material without resorting to expensive trial-and-error is a limitation for its extensive industrial use.Understanding the nature of the interactions between the layered material and the solvent would allow for the preemptive screening of effective solvents for the exfoliation of layered materials.As well as aiming to maximize monolayer yield, there are other restrictions on the choice of solvent.For example, many of the better solvents for graphite exfoliation, including NMP, DMA and DMF, are toxic and facing restrictions by the European Chemicals Agency, rendering them unsustainable for future industrial use [22].Devising a screening tool to identify non-toxic, low-boiling point solvents which do not comprise on performance would be extremely useful [23,24].Ideally these screening tools would be generalizable to other two-dimensional materials beyond graphene.
One of the first attempts to screen for effective solvents found that the concentration of dispersed graphene was maximized for those solvents with a surface energy or tension similar to that of graphite, following the traditional surface wettability argument [12].Although this method has had some successes, there are also examples where it breaks down [20].Hansen solubility parameters, characterized by three intermolecular interactions, namely hydrogen bonding, polar and nonpolar (dispersive), have also been used as screening parameters.Graphene was found to have Hansen dispersion, hydrogen bonding and polar solubility parameters of δ D = 18 MPa 1/2 , δ P = 9.3 MPa 1/2 , and δ H = 7.7 MPa 1/2 , respectively [20].Searching for solvents with solubility parameters close to those of graphene led to the prediction and confirmation of cyclopenanone and cyclohexanone as excellent solvents of graphene.However, here again, there are some examples where this screening technique fails -dimethyl phthalate has very similar Hansen parameters to cyclopentanone, yet is a poor graphene solvent [20].A peculiarity of the Hansen solubility parameters attributed to graphene is the non-zero polar component (δ P = 9.3 MPa 1/2 ).The mechanism which could lead to this is so-far unknown.Contribution from edge-sites or the unintentional chemical functionalisation of the basal plane were touted as a potential physical mechanisms, yet there is little experimental evidence for either [25,26].All of the above suggests that empirical solubility models, which consider only macroscopic solution thermodynamics, are blunt tools, missing some information about the stabilization mechanism [27,28].These methods assume that kinetic effects driving reaggregation play a minor role.There is growing acceptance that this is unlikely to be true and that explicit structural and electronic interactions between the solvent molecules and the solute are important in the stabilization of the exfoliated monolayers [29][30][31][32].
Here, we systematically study the origin of the solvent dependence of the stabilization of exfoliated monolayers.We address several open questions: Is it sufficient to consider only energetic effects via the enthalpy of mixing when choosing screening descriptors, what is the physical mechanism leading to the experimental finding of a non-zero polar component of the Hansen solubility parameter of graphene and, and how important are kinetic effects, such as reaggregation, to the stabilization of graphene in a particular solvent?

Choice of solvents
The twelve solvents included in this investigation are listed in Table 1.They were chosen to include those reported by Hernandez et al. [33] to result in a wide range of graphene concentrations after exfoliation and include cyclopentanone (best performing solvent, polar aprotic), N-methyl-pyrrolidone (commonly used effective solvent, polar aprotic), benzaldehyde (poor performance, polar aprotic), bromobenzene (intermediate performance, slightly non-polar) and toluene (poor performance, non-polar).All others have intermediate exfoliation capabilities or polar properties.Note that while the concentration of exfoliated graphene in individual solvents has been increased considerably since the initial work of Hernandez et al, it is, to the best of our knowledge, the only available comparative study of the effectiveness of a wide range of solvents for graphite exfoliation.

Molecular Dynamics
Molecular dynamics calculations are performed using the Large-scale Atomic / Molecular Massively Parallel Simulator (LAMMPS) [34].As the interaction between graphene and isolated solvent molecules was found to be primarily attractive dispersion forces [35,36], it can be modelled accurately using classical van der Waals (vdW) force fields.The graphene monolayer is modeled as uncharged vdW spheres and are described by interaction parameters originally reported by Steele et al. [37] and later used to describe the interaction between graphene and molecular adsorbants [30,[38][39][40].The molecules are described by the All Atom Optimized Potentials for Liquid Simulations (OPLS-AA) potentials, which were obtained from LigPargen web server [41][42][43].The pair interaction coefficients between graphene and the molecules are obtained using Lorentz-Berthelot rules [44,45].Binding energies calculated using these potentials agree with experiment for small organic molecules adsorbed on graphene [35].The initial solvent structure was created using packmol [46] in a box size of 105.91 Å × 106.65 Å × 120.0 Å to reproduce the experimental roomtemperature densities (see Table B1).The graphene monolayer was then included in the simulation box by increasing the box size in the direction normal to the plane by (2 × 3.5) Å.The python package pymatgen was used to generate the LAMMPS structural data file [47].The systems were annealed from 600 K to 300 K for 50 fs with an integration timestep of 0.01 fs and then equilibriated at 300 K for 1 ns with an integration timestep of 1 fs.The system was further simulated with the NVE ensemble for 2 ns, of which 2000 samples from the last 1 ns were used to generate statistical distribution functions.For the duration of the complete simulation the graphene layer is held fixed.

Thermodynamics of Mixing
The Helmholtz free energy of mixing, ∆A mix , is an indication of the favorable formation of a solvation structure.As ∆A mix = ∆H mix − T ∆S mix , where ∆H mix and ∆S mix are the enthalpy and entropy of mixing, respectively, a requirement of mixing is that ∆H mix < T ∆S mix .Note that the PV term of the enthalpy is negligible for liquids and solids so it can be neglected.In this case, the internal energy of mixing at 0 K will be equal to the enthalpy of mixing at 0 K. ∆A mix was determined from molecular dynamics simulations using achemical free energy methods [48].Such methods simulate a series of nonphysical intermediates to calculate the free energy of transferring a solute (here graphene) from vacuum (state 0) to solution (state 1).A continuous variable λ parameterizes the path between state 0 and state 1. ∆A mix is calculated here using the finite difference thermodynamic integration (FDTI) method [49][50][51].MD simulations are carried out for a discrete set of λ values and the free energy of mixing is calculated as: where ⟨. ..⟩ indicates an ensemble average, w i are integral weights and δ is the finite difference parameter which satisfies (λ i − λ i−1 ) >> δ.To avoid numerical instabilities the Leonard-Jones potential is modified to be a function of λ.The system is simulated for 10 ns in the NVT ensemble at 300 K.The parameter λ is incremented in 40 steps (every 250 ps) with the finite difference parameter δ = 0.001.For each λ, the first 100 ps are discarded for equilibration and the ensemble average is calculated by sampling every 20 fs from the last 150 ps.The enthalpy and entropy contributions to the total free energy are extracted from the slope and yintercept of ∆A/T versus 1/T , respectively, where the free energy is calculated at 290 K, 300 K and 310 K.

Diffusion Coefficients
The diffusion coefficient of solvent molecules confined in a region [a, b] parallel to the graphene sheet, D xx , was determined from MD trajectories using the method described in Liu et al. [52] and Agosta et al. [53].D xx is calculated as: where the terms on the right hand side can all be calculated using standard molecular dynamics calculations.⟨∆x(τ ) 2 ⟩ {a,b} is the mean square displacement of particles in the direction parallel to the surface in a region {a, b} in a time window τ and P (τ ) is the probability that a particle remains in the region {a, b} during that time, i.e, the survival probability.⟨∆x(τ ) 2 ⟩ {a,b} is extracted from MD trajectories via: where ζ(t, t + τ ) is the set of all particles that stay in the volume a < z < b from time t to t + τ , N (t) is the number of particles in the region {a, b} at time t, and T is the number of time steps to average.The survival probability is calculated as: where N (t, t + τ ) is the number of particles that stay within the region {a, b} between times t and t + τ .Note that the parameters in the calculation of diffusion coefficient are τ , T , and the definition of the region between a and b.Here, the equilibrated structure was simulated with the NVE ensemble for T = 150 ps and the center of mass of each molecule recorded at each time step.A time window (τ ) of 35 ps is used to calculate the diffusion coefficient.a and b are chosen to encompass the entire first solvation shell of each system.

Results
In Section 3.1, we show that the free energy of mixing is always negative regardless of the solvent due to a large enthalphic contribution for solvating graphene.A significant entropic penalty, however, suggests that the solvent must undergo significant ordering in the presence of graphene.Section 3.2 discusses this structural and orientational ordering.Finally, Section 3.3 presents the diffusion coefficients of solvent molecules in the first solvation shell and discusses how these molecules play a major role in preventing the reaggregation of graphene sheets.

Free Energy of Mixing
It was previously shown, using ab initio calculations, that the interaction between an isolated gas phase solvent molecule and a pristine graphene monolayer is van der Waals in nature, with no charge transfer involved.This is independent of the type of molecule, and the binding strength of a solvent molecule on graphene is not correlated with its ability to exfoliate graphite [36].To determine if instead the collective behaviour of the molecules near the surface of a graphene monolayer is the determining quantity, we calculate the Helmholtz free energy of mixing, ∆A mix , of graphene in a solvent.If ∆A mix < 0, then the mixture is thermodynamically more stable than each of the two components separately.The Helmholtz free energies of mixing for graphene at room temperature are shown in the first column of Table 1 for a variety of solvents.At 300 K, ∆A mix is always negative, irrespective of the solvent, suggesting the favorable mixing of monolayer graphene in each of the solvents.
It is most favourable in quinoline (∆A mix = −11.08MJ/mol) followed by bromobenzene (∆A mix = −11.41MJ/mol).It is least favourable in cyclohexanone, at −7.57MJ/mol.The magnitudes of ∆A mix reported here are consistent with those reported by Oyer et al. for the free energy of mixing of graphene in benzene and hexaflourene [54] and by Mukhopadhyay et al. for the free energy of mixing of graphene in water and DMSO [55].
The decomposition of the Helmholtz energy into the enthalpic and entropic contributions are given in the second and third column of Table the other solvents at −9.57MJ/mol, while the entropy contribution is positive at 0.83 MJ/mol.However, this unusual behavior is not evident in the overall Helmoltz free energy of mixing which is similar to that of the other considered solvents, at −10.41 MJ/mol.
In general, no correlation can be found between the magnitude of ∆A mix and the experimentallydetermined graphene concentrations after exfoliation.
As such, the Helmoltz free energy of mixing cannot be used as a mechanism of screening for more effective solvents.Instead, the calculated free energies indicate that the solvent-dependent exfoliation is not energetically bound and other dynamical effects should be considered.

Solvation Shell Formation
The entropic penalty suggests that there must be some solvent reconstruction or ordering that occurs due to the presence of graphene.To understand this, the location and orientation of solvent molecules close to the graphene layer are extracted from molecular dynamics calculations.The position distribution function g(z) and the angle distribution function α(z) are shown in Fig. 1.The position distribution function g(z) is given relative to the bulk solvent number density so that a value greater than 1 indicates an accumulation of molecules relative to the bulk solvent, and a value less than 1 indicates a depletion of molecules relative to the bulk.α(z) is the angle made by the normal to the plane of the molecule with respect to the z axis, where the plane is defined as containing at least 3 atoms of the molecular ring system.
Irrespective of the polar nature of the molecule, distinct solvation shells are formed next to the graphene layer as a result of surface confinement, molecule -molecule interactions and molecule -graphene interactions.The first solvation shell (grey shaded area) can be recognized as a peak in g(z) followed by a deep trough near the surface of graphene.The peak associated with the first solvation shell is also the sharpest, indicating that the highest degree of transverse ordering occurs for those closest to graphene.The first solvation shell has a complex structure in most cases due to the adsorption geometry of the molecules on the graphene surface.This is particularly evident for cyclopentanone (Fig 1(a)) and NMP (Fig 1(f)) where the first solvation shell is composed of two distinct peaks.The reason for this can be found by looking at α(z): the solvent adsorbs on graphene with two distinct angular orientations -one with the molecular plane lying approximately parallel to the surface (α = 0 • ) and one with the plane almost perpendicular (α = 90 • ).As the distribution function is based on the centre of mass of the molecules, this manifests as a double peak structure in g(z).As a result, the dipole moments of the molecules also show two distinct polar angles (Θ(z)), whereas the azimuthal angle (ϕ(z)) shows no distinct orientation preference indicating no in-plane ordering of the dipoles (see Fig A1).This behavior was previously found in molecular dynamic calculations involving NMP interacting with both graphene and carbon nanotubes [39,56,57].In other cases, the distinction between the two peaks is not as clear.For example, the first solvation shell of cyclohexanone (Fig 1(b)) has a small peak corresponding to those molecules located closest to graphene and orientated parallel to it, but the remainder of the molecules in the first solvation shell do not exhibit an angular preference.This is also the case for NFP, bromobenzene, chlorobenzene, dioxolane, quinoline, acetone and toluene.For the case of benzonitrile and benzaldehyde the first solvation shell is comprised almost entirely of molecules lying parallel to the graphene layer, with no second peak visible.There is some experimental evidence for the confinement of solvent molecules near the surface of graphene -Arunachalam et al.
showed that the NMP molecules near the surface of graphene show a reduction in rotational degrees of freedom using rotating frame Overhauser effect spectroscopynuclear magnetic resonance technique [58].
To visualize this orientational ordering, Fig 2 shows a snapshot from the molecular dynamics simulation showing only the first solvation shell of benzonitrile and NMP.All molecules belonging to the very first peak are shaded in red, and all others are shaded in green.As expected in the case of benzonitrile, it can be seen that essentially all molecules in the first peak are orientated parallel to the plane of graphene (α ≈ 0 • ) with very few molecules adsorbed perpendicularly.In contrast, a small but significant number of NMP molecules are adsorbed perpendicular to the sheet.These molecules appear to fill all the gaps that exist in between the flat-lying molecules.To determine if any long-range lateral order is present, we calculated the pair correlation functions for the geometric centers for the molecules in the first solvation shell.These are shown in Fig 3.In all cases, there is a local ordering around the center (0,0) due to the presence of an exclusion zone around the molecule.This ordering then decays rapidly moving away from the center.The strongest ordering is seen in case of cyclopentanone, cyclohexanone and 1,3-dioxolane.Intermediate ordering behaviour is seen in case of NFP, benzonitrile, NMP, quinoline, acetone and toluene, with no ordering beyond the first shell seen in case of bromobenzene, chlorobenzene and benzaldehyde.This is in agreement with the work of Terrones et al. who found that while a monolayer of NMP molecules confined between two graphene sheets exhibit longrange hexagonal ordering, this ordering is lost as the thickness of the molecular layer increases [57].
In response to the formation of the first molecular layer, the rest of the solvent then reorganizes to form further Table 2. Data extracted from the pair distribution function g(z) shown in Fig. 1. h min is the distance between graphene and the solvent atoms closest to it.N shells indicates the number of solvation shells present where the final solvation shell is defined to that peak which deviates from the bulk value by at most 5%.Depth is the distance of last solvation shell from graphene.Period indicates the average periodicity of the solvation shells.
solvation shells, as seen in Fig 1 .As expected, the amplitude of the solvation shell peaks decay as the influence of the solute wanes, finally approaching a constant value of 1, indicating bulk solvent behaviour.
The distance between graphene and the nearest solvent atoms (h min ), the number of solvation shells (N shells ), the depth to which the solvation shells extend away from graphene, and the average distance between two consecutive solvation shells are given in Table 2.
As the distance between graphene and the atom of the molecule nearest it is determined by the vdW interaction, it is very similar across all solvents regardless of the molecular orientations, varying from 1.99 Å for bromobenzene to 2.07 Å for benzaldehyde.
The number of solvation shells varies from a minimum of 3 in the cases of bromobenzene, benzonitrile, chlorobenzene and benzaldehyde to a maximum of 6 for cyclohexanone and 1,3-dioxolane, where in counting the total number of solvation shells we have defined the last shell as the final peak that deviates from bulk by at most 5%.The average distance between adjacent solvation shells ranges from 4.25 Å for benzonitrile to 5.25 Å for quinoline.
Despite these differences, it is clear from Fig 1 that the solvation structure cannot be used to explain the observed phenomenon of solvent-dependent graphene stabilization.For example, there is no appreciable difference between the solvation structure of NMP and toluene, one of the best and worst solvents for the exfoliation of graphite, respectively.However, the existence of these solvation layers can explain the origin of the non-polar Hansen solubility parameter attributed to graphene.Rather than to graphene alone, the solubility parameters can be attributed to an effective solute comprised of graphene and its first few solvation shells.Due to partial charges on the atoms of the first solvation shell, the effective solute appears to have surface charges, resulting in a designation of a non-zero polar solubility parameter.

Molecular Diffusion
In the analysis up until now, we have assumed that the process of stabilization is thermostatically driven, and kinetic effects such as sedimentation play a minor role.This may not be necessarily true, and in certain conditions may be the most important consideration.
If the solvent-dependence of graphene stability in a dispersion is determined by reaggregation effects, then screening descriptors based on energetic considerations alone will not be not appropriate.Long-term stability can only be obtained if the dispersion is kinetically stabilized against reaggregation when layer collision occurs.If the layers are far apart, this will be influenced by the viscosity of the solvent [24,59].
When the layers are close together, the diffusion coefficient of molecules confined close to the graphene sheet will determine the ease at they can be ejected to facilitate reaggregation.Indeed, the addition of surfactants or polymers to graphene dispersions to prevent reaggregation operate on this principle [60,61].
To determine the extent to which solvent molecules can kinetically block reaggregation, we calculate the diffusion coefficient in the direction parallel to the graphene surface of the molecules in the first solvation shell.
The results are shown in Fig. 4 as a function of the experimental graphene concentrations taken from Ref. [33].The highest parallel diffusion coefficient is found for acetone, at 0.35 Å 2 /ps.This is followed by toluene (0.18 Å 2 /ps), dioxolane (0.16 Å 2 /ps) and chlorobenzene (0.16 Å 2 /ps).It is notable that all of these are poor solvents for the exfoliation of graphite, with the worst solvent corresponding to that with the largest parallel diffusion coefficient.As the diffusion constant decreases, the experimental concentration increases, with the best solvent (cyclopentanone) having a parallel diffusion constant of 0.085 Å 2 /ps.However, in the very-low diffusion regime the concentration does not follow this upward trend.As the molecular hindrance to re-aggregation is higher in this regime, one would expect that the exfoliated layers would remain in the dispersion resulting in an increased monolayer concentration.Instead, the concentration decreases again.
While this requires further investigation, we hypothesis that the concentration of graphene is now dependent on the ability of that solvent to faciliate the exfoliation process itself, as opposed to the stabiliziation process.When graphene is dispersed in a solvent with a low surface diffusion coefficient, such as benzaldehyde, the reaggregation rate will be low and the concentration of graphene layers will be maintained at its initial concentration.
In contrast, when the diffusion coefficient is high, the graphene concentration will always be low due to a high rate of reaggregation.Such solvents will result in a low concentration dispersion even if a large number of monolayers are initially exfoliated due to efficient reaggregation.As a consequence then, solvents with low surface diffusion coefficients, such as quinoline and benzaldehyde, are good candidate for the optimization of the initial separation process.If the initial concentration of graphene in these solvents can be increased, the low diffusion rate in the first solvation shell means that reaggregation will be hindered and the high concentration maintained.

Conclusion
We find that there is a molecular-level structural and orientational ordering of solvent molecules at the graphene-solvent interface which extends up to 30 Å away from the graphene layer.These solvation layers form in all cases and behaves similarly irrespective of the solvent polarity.
We speculate that the non-polar Hansen parameter obtained by fitting the experimental data can be attributed to the formation of the first solvation shell as the molecules in the first solvation shell and the graphene form an effective solute.The formation of these solvation layers is spontaneous as indicated by the free energy and its decomposition in enthalpy and entropy contributions.This independence of the solvation structure on the nature of the solvent molecule suggests that non energetic contributions are responsible for the observed solvent-dependence of graphene stabilization, and specifically the kinetics of reaggregation.We propose that the diffusion coefficient of molecules in the first solvation shell is an appropriate property to determine the probability of reaggregation in a solvent.We note that in the low-diffusion limit, this appears to break-down, at least for the solvents investigated here.This origin of this requires further work, but could be related to the role of the solvent in the initial separation step.We predict that the concentration of graphene exfoliated in solvents with low diffusion coefficients parallel to the graphene plane, such as benzaldehyde and quinoline, can be enhanced dramatically by careful tuning of the experimental sonication parameters, as any produced monolayers will be prevented from reaggregation by steric effects due to the slowly diffusing solvent molecules on the monolayer surface.

Figure 1 .
Figure 1.Pair distribution functions.(Top panel) g(z) is the position distribution of the center of mass of the solvent molecules.g(z) is given relative to the bulk solvent number density.(Bottom panel) α(z) is the probability distribution function of the angle made by the normal to the plane of molecule with the axis perpendicular to the plane of graphene.The graphene monolayer is located at z = 0 Å.

Figure 2 .
Figure 2. A snapshot of molecules in first solvation shell of benzonitrile (left) and NMP (right) taken from the molecular dynamics simulation.Red shading indicates molecules belonging to the first peak in the solvation shell and green shading indicates the rest of the molecules.

Figure 3 .
Figure 3. Pair correlation function for the center of geometry of the solvent molecules present in the first solvation shell.

Figure 4 .
Figure 4.The diffusion coefficient of solvent molecules in the first solvation shell in the direction parallel to the surface of graphene as a function of the experimental concentration of graphene C G (µg/ml) in that solvent.Experimental data is taken from Herndandez et al. [33].

Figure A1 .
Figure A1.Pair distribution functions.(Top panel) g(z) a histogram of the center of mass of the solvent molecules normalised by ( N L dz) where N is the number of molecules in the unit cell, L is the length of the unit cell and dz is the bin width.The graphene monolayer is located at z = 0 Å. (This is the same as presented in Fig. 1).θ(z) is the probability distribution function of the polar angle of the molecular dipole, where 90 • corresponding to the molecular dipole in the plane parallel to graphene and 0 • (180 • ) corresponds to the molecular dipole in the plane perpendicular to the plane of graphene.(Bottom panel) ϕ(z) is the probability distribution function of the azimuthal angle of the molecular dipole, indicating the orientation of the dipole in the plane parallel to the plane of graphene, with 0 • (180 • ) corresponding to the orientation along positive (negative) x-axis, 90 • (−90 • ) corresponding to the orientation along the positive (negative) y-axis.

Table 1 .
1, respectively.In almost all cases there is a significant entropic penalty for solvating graphene, with T ∆S mix lying between −17.13 MJ/mol and −12.07 MJ/mol.This entropic penalty is compensated by a large enthalpy change ranging between −28.32 MJ/mol in quinoline to −21.10 MJ/mol in toluene, driven by strong van der Waals interactions.Bromobenzene is an exception to this: ∆H mix is much smaller in magnitude than for Solvent ∆A mix (MJ/mol) ∆H mix (MJ/mol) T ∆S mix (MJ/mol) The Helmholtz free energy of mixing, ∆A mix , the corresponding enthalpy, ∆H mix and the entropy term calculated at room temperature (300 K), T ∆S mix .These values are determined for a graphene area of 112.953 nm 2 .