Gas–liquid phase equilibrium of a model Langmuir monolayer captured by a multiscale approach†

The paper exemplifies a multiscale approach towards the determination of effectively two-dimensional gas–liquid phase transitions of Langmuir monolayers. Based on information gathered from molecular dynamics simulations of three-dimensional model Langmuir chain-like surfactant monolayers we calculate an effective lateral interaction potential between projected surfactant center-of-masses (red coins) using a force-matching scheme. The obtained coarse-grained interactions enter Gibbs ensemble Monte Carlo simulations of the corresponding two-dimensional system, and are further used to derive implications of two-dimensional density functional theory based on Weeks–Chandler–Anderson perturbation theory. As featured in:


Introduction
Langmuir monolayers are surfactant films formed by amphiphilic molecules such as lipids at an air-water interface.2][3][4] There has been a special interest in Langmuir monolayers of fatty acids and phospholipids as model systems for lung surfactants 5 and cell membranes which can be considered as weakly-interacting monolayers. 6epending on surfactant architecture and solvent type, 4,7 monolayers have a very different but rich phase diagram which is mainly characterized by the density of surfactants, chain alignment and spatial ordering 1 at the interface.At extremely low surface concentrations, they exist in a 2D gaseous (G) phase.
,10 Langmuir monolayers have been mainly characterized experimentally by their 2D surface pressure-area (P-A) isotherms.Different regions in the P-A isotherms are attributed to various phases and their coexistence regions of the corresponding Langmuir monolayers: a horizontal line in the P-A isotherm is interpreted as the first-order G-LE phase coexistence regime, 9,11 or a G-LC phase transition for temperatures below the G-LE-LC triple point, 12,13 whereas a non-horizontal line is found to exist for coexistence of condensed phases. 7The formation of 2D surface micelles, 2,3,7,14,15 which has been observed experimentally by using e.g.fluorescence microscopy, is found to be the main reason for the non-horizontal line in P-A isotherms. 1,16][18][19] The determination of precise phase transition boundaries and a critical point (if present) of the Langmuir monolayers is extremely difficult if not impossible from experiments. 2,9,20heoretical frameworks such as Landau theory of phase transitions have been used with some success in predicting condensed phases and coexistence boundaries of Langmuir monolayers 21 qualitatively.Molecular dynamics (MD) and Monte Carlo (MC) simulations 1,2,8,10,[22][23][24] of the full 3D system including the liquid phase do not only provide more quantitative information but have also broadened our insight into the structure and phase behavior of Langmuir monolayers.Current state-of-the-art atomistic simulations 25 were so far successful in predicting the structure of gas and solid phases of phospholipid Langmuir monolayers (76 nm 2 in size, with 76 to 144 surfactants at each interface).However, when it comes to the direct simulation and determination of phase boundaries and the critical point of such complex systems, they do not provide us with sufficiently reliable data since their finite sizes affect the phase segregation at the interface dramatically. 9,22specially for the G-LE phase transition taking place at a very low surface concentration one needs to consider very large 3D systems with large enough interfacial area to accommodate a sufficient amount of surfactants to avoid finite-size effects. 26his issue can be addressed by the Gibbs-ensemble Monte Carlo (GEMC) technique 27 which can be used to determine liquid-vapor phase equilibria even by considering a relatively small number (E1000) of particles.This technique has been applied successfully to simple liquids such as 2D 28 and 3D 27 Lennard-Jones fluids.
To the best of our knowledge, there has been so far only a handful of GEMC studies of Langmuir monolayers in the literature. 9,29Despite the fact that solvent conditions can completely change the phase behavior of the monolayer, 4 solvents were treated only implicitly due to technical difficulties, and a configuration-biased Monte Carlo (CBMC) method had to be used to move bead-spring chain surfactants with a reasonable acceptance rate. 9As an alternative approach we have recently developed a multiscale coarse-graining (MSCG) approach which maps the entire 3D system to a 2D system consisting of surfactant center-of-masses adsorbed at the fluid interface. 26ere, we apply this approach to a 3D model reference system 8 consisting of linear bead-spring chain surfactants spread at a vapor-liquid interface, and calculate effective pairwise interaction potentials at relatively small surface concentrations.The effects of solvents and the chain structure are taken care of in this coarse-grained (CG) interaction potential. 26For the corresponding 2D systems of surfactant center-of-masses described by this effective interaction potential, similar to a simple 2D fluid, we use density functional theory 30 (DFT) and GEMC simulations 31 to determine G-LE phase boundaries along with thermodynamic properties such as pressure and chemical potential at the coexistence region.By fitting the results of the GEMC simulations to the critical behavior of the Ising model 32 (with a different critical exponent b, however) we estimate the critical temperature, and by using the law of rectilinear diameters, we estimate a value for the critical density. 31he remainder of this paper is organized as follows.In the next section, we discuss the details of the reference MD simulations and present selected simulation results.We then perform multiscale coarse-graining of the reference systems and extract effective interaction potentials.We verify the accuracy of our force matching scheme by comparing radial distribution functions obtained via independent MC simulations of the 2D CG system with those of the reference simulations.Next, we use the results of the CG interaction potentials to perform DFT calculations and GEMC simulations and determine 2D G-LE phase boundaries within the monolayer.Areal densities, chemical potentials and surface pressures of the G and LE coexisting phases are calculated from both theory and simulation and compared with each other.We further use DFT calculations to determine equilibrium density profiles across the G-LE dividing line and corresponding line tensions.Finally, we summarize our results in the Conclusions section.Supporting details on multiscale coarse-graining, MC simulations and DFT calculations are discussed in the Sections S2 and S3 (ESI †).

Reference simulations
We consider a simple well-known molecular model which has been introduced and investigated for its phase behavior by Tomassone et al. 8 The 3D reference system consists of monoatomic solvent particles, namely w (water-like), and linear surfactant molecules HT 5 with one hydrophilic head (namely H) and five hydrophobic tail beads (namely T).All pairs of particles (also denoted as ''beads'') interact via a truncated Lennard-Jones (LJ) potential, shifted by V cut ij so that V ij (r cut ) = 0. Coefficients s and e denote the preferred diameter of particles and the corresponding potential well depth, while i, j A {H, T, w}.All beads are considered to have the same mass m and the same preferred size s.The cut-off distance, r cut , is set to 2.5s for all pairs in the system.The coefficients c ij control the attractive branch of the LJ potential and hence the affinity between a pair of particles of type i and j.The corresponding numerical values c ij for different particle pairs are given in Table 1.These coefficients are tuned to model the insoluble case where surfactants form a Langmuir monolayer by residing preferably at the vapor-liquid interface with their head attached to the liquid sub-phase and their tail pointing out toward the vapor phase. 8We ensure surfactant chain connectivity by considering the additional finitely extendable nonlinear elastic (FENE) bond potential between successive beads as 8 where R 0 = 1.5s is the maximum bond length and k = 30e/s 2 is the spring constant.The temperature of the 3D reference system was fixed at T = 0.9e/k B within the vapor-liquid coexistence region of the solvent phase.The full bond potential V FENE (r) + V LJ (r) has a minimum at r E 0.961s which is close to the mean equilibrium bond length r E 0.969s resulting from the Boltzmann distribution at the chosen temperature.
We have used the open-source molecular dynamics simulation package, LAMMPS, 33 to perform 3D reference simulations in the NVT ensemble.A periodic simulation box of lengths L x = L y = 50.5sand L z = 150.5s,containing a total number of 9 Â 10 4 particles of type w, H and T has been used for all simulations.An integration time step of 0:007s ffiffiffiffiffiffiffiffi m=e p was used to integrate modified Hamilton's equations of motion as the temperature was fixed by using a Nose ´-Hoover 34 thermostat with a damping time of 0:7s ffiffiffiffiffiffiffiffi m=e p .A smaller system with the same interaction potentials at the same state point at T = 0.9e/k B has been studied previously. 8o prepare a Langmuir monolayer, we generate first a bare vapor-liquid interface of pure w particles and then place HT 5 surfactants at the interfacial region.To this end, we first simulate a pure liquid slab of w particles at the state point T = 0.9e/k B , r E 0.85s À3 in a periodic simulation box with the same L x and L y and a smaller dimension of 41.5s in the z-direction.After 5 Â 10 4 time steps, this liquid slab is translated into the middle of the main simulation box (with L z = 150.5s)and is relaxed for additional 10 5 time steps until reaching equilibrium with two vapor-liquid 2D interfaces normal to the z-direction.In order to introduce surfactant molecules into the system, six adjacent solvent particles (within a neighborhood of 1.4s radius) are selected randomly from either of the interfacial regions.These randomly selected w particles are then connected together with FENE bonds and their identities are changed into H and T, such that the total number of particles remains unchanged after a surfactant insertion.We continue this procedure until placing a number of 2N int surfactant molecules at the two interfacial regions and hence reaching a desired surface concentration G = N int /L x L y .The MD simulation is then continued for another 1.5 Â 10 5 time steps and thermodynamic properties such as pressure and temperature are monitored to make sure the system has equilibrated before the sampling interval.We continue the simulation for another 1.5 Â 10 5 time steps and collect samples for statistical analyses once every 30 time steps.
A simulation snapshot together with the density profiles in the system at surface concentration G = 0.118/s 2 are shown in Fig. 1.Surfactants form Langmuir monolayers normal to the z-direction at the vapor-liquid interfaces (Fig. 1a) with their heads (yellow beads) immersed in the liquid sub-phase and their tails (red beads) pointing towards the vapor phase.This is more clearly visible from Fig. 2 where we provide top views of  the system at different surface concentrations.The surfactants are seen to form rather homogeneous monolayers with a weak tendency to form clusters. 8 These qualitative observations can be quantified by number density profiles for different species as shown in Fig. 1b, where n s (red dashed-line) represents the number density of the surfactant center-of-masses.The surfactant center-of-mass profile is peaked at the interfacial regions and vanishes in the two bulk phases, verifying the observation that the system is forming insoluble Langmuir monolayers. 8The water density n w profile indicates uniform densities in the bulk phases (with average values of r l = 0.8569 AE 0.0002s À3 and r v = 0.0063 AE 0.0002s À3 ) and is slightly peaked at both interfaces due to the higher value of the attraction coefficient c wH than c ww beads.The water particles thus prefer H particles as neighbors rather than their own species.This fact can be inferred from radial distribution functions of H-w and w-w pairs and hydration numbers of H groups and w particles given in Fig. S1 of the Section S1.1 (ESI †).The peak in the total particle number density n p is even more pronounced due to the adsorption of the H and T monomers at the interface plane.
Surface tension g is calculated employing its standard definition. 32For a 3D periodic system with two planar interfaces normal to the z-direction, the expression takes the form where hÁ Á Ái denotes the canonical ensemble average and P xx , P yy and P zz are diagonal components of the pressure tensor.The surface tension isotherm of the system evaluated from eqn (3) via its virial expression for constant T and varying surfactant concentration is depicted in Fig. 3.The data shows that g decreases by increasing the surface concentration in the range we have examined.We note that the values of the surface tension are slightly different from what was reported by Tomassone et al. 8 for a system with the same interaction potentials at the same state point.The discrepancy is most probably due to finitesize effects, the number of particles in our present study is almost one order of magnitude larger.Furthermore, in contrast to these authors we do not observe a statistically relevant plateau in the surface tension-area isotherm whose presence would signal the phase coexistence regime.The observation of a phase transition is known to be impossible in the simulation of smallsized monolayers, 22 especially in full 3D simulations where one needs to consider a relatively large number of solvent particles to just create a gas-liquid or liquid-liquid interface.
Next, we analyze the surfactant monolayer itself more closely by investigating the surfactant center-of-mass density profile and bond uniaxial orientational order parameter S A [À0. 5, 1] in the vicinity of the 2D interface, where the latter is defined as Here P 2 is the second Legendre polynomial and y the polar angle between an individual bond within a surfactant molecule and the z-direction normal to the 2D interface.The density profiles in Fig. 4a show that the surfactant center-of-masses are confined within a narrow region at the interface with quite a small thickness roughly of the size of the stretched surfactant end-to-end distance (E5s).Furthermore, the bond order parameter S, illustrated in Fig. 4b, shows a weak orientational ordering within the surface concentration range we have examined.Therefore, in this regime we neglect orientational ordering effects and apply our recently developed approach 26 to map the entire 3D system to an effective 2D simple fluid of surfactant center-of-masses.
3 Coarse-grained interaction potential The 3D reference simulations can be used to extract the effective interaction between surfactants within the monolayer constituting the 2D interface.To this end we first consider a simple free energy argument that involves only an effective hard disc diameter for the adsorbed surfactants at the interface plane.Afterwards, we discuss in detail the multiscale coarsegraining procedure to systematically map the whole 3D system to an effective 2D system described by an effective interaction potential.
Fig. 3 Surface tension as a function of surface concentration for the 3D reference system.Dash-dotted line is just a guide for the eye.We use the surface pressure-area isotherm to determine an effective disc diameter d of surfactants adsorbed at the 2D interface.Fig. 5 (symbols) shows the surface pressure P as a function of area per molecule, S = A/N int = 1/G, defined by where g 0 denotes the surface tension of the bare interface and N int is the total number of surfactants adsorbed at the assumed flat interface with area A = L x L y .In order to rationalize Fig. 5, we employ a simple argument first advocated by Pieranski 35 for colloids and later used for polymers. 36,37Each effective surfactant molecule removes an area A 1 = pd 2 /4 from the bare vapor-liquid interface, where d is an effective diameter.The resulting change of free energy can be approximated by and a corresponding surface tension thus approximately given by The surface pressure resulting from eqn ( 5) and ( 7) is found to be inversely proportional to the area per molecule, Despite its simplicity, eqn (8) can describe the results of the simulated surface pressure-area isotherm with an effective diameter d = 1.41 AE 0.09s reasonably well, as shown by the solid line in Fig. 5.We note that this simple free energy analysis neglects (possible) attractive interaction potentials between the hard discs.In order to consider that, we fit the surface pressurearea data to the commonly used Volmer equation of state 17 for 2D gaseous insoluble monolayers as where P* is the cohesion pressure which accounts for the attraction between hard discs.As it is shown in Fig. 5 by the dash-dotted line, eqn (9) describes the simulation results very accurately with an effective diameter of d Vol.= 2.2 AE 0.1s and a cohesion pressure of P* = 0.002 AE 0.005e/s 2 as fitting parameters.This fitting estimates a somewhat larger hard disc diameter and suggests the existence of an attractive interaction potential between the hard discs.We will compare the values of d and d Vol.later to an alternative definition of an effective hard core diameter.

Multiscale coarse-graining
We now discuss the multiscale coarse-graining of the 3D system containing the monolayer and its surrounding to an effective 2D system carrying the potentially phase-separating monolayer.
More details are available in our previous study. 26Here we focus on aspects of relevance for the present manuscript.Within this approach one defines N int CG sites R I as surfactant center-of-masses projected onto the xy-plane described by a mapping M R I , more formally, where r n denotes the position vectors of all particles in the 3D reference system and N int is the total number of surfactants adsorbed at one monolayer.Since the mass of all beads within the surfactants was chosen to be identical, the mapping operator takes the form where (ı ˆ, j ˆ, k ˆ) denote Cartesian unit vectors with k ˆthe surface normal, r s Ij denotes the position vector of the jth bead within the Ith surfactant and N b is the total number of beads of the surfactant chain.Note that the mapping (11) ignores the height of the surfactant center-of-mass which is justified in view of the small interface width (see Fig. 4a).We follow ref.38 and require consistency between coarse-grained and molecular-level descriptions of the corresponding configurations, R Nint = {R 1 ,. .., R Nint } and r n in the canonical ensemble, where U(R Nint ) and u(r n ) are the CG and reference total potential energy, respectively, and b = 1/k B T. eqn (12) implies that (i) the coarse-grained potential U(R Nint ) can be determined up to an additive constant by the reference potential and the mapping operator, and (ii) it is a function of temperature and surface concentration through its dependence on N int .We follow common practice 39 and assume a pairwise interaction potential U CG , which due to symmetry within the interface plane is a central potential resulting in a radial pairwise force strength as between CG sites.The total potential energy of the coarse-grained system is therefore given by and the coarse-grained force exerted on the Ith projected surfactant center-of-mass reads The force matching scheme 38,39 amounts to determining the coarse grained forcefield F CG (r) from underlying 3D reference simulations by minimizing the residual where hÁ Á Ái denotes the canonical expectation value and f I is the net force exerted on the Ith surfactant center-of-mass, projected onto the xy-plane.Implementation details of the force matching scheme are provided in the Section S1.2 (ESI †).Schematic Fig. 6 shows the entire multiscale coarse-graining process: the mapping operator takes the HT 5 surfactants and projects their center-of-masses onto the xy-plane; the corresponding forcefield is subsequently calculated through the force matching scheme.The CG interaction potential obtained from integrating F CG (r) is shifted such that U CG (r) vanishes at r = 10s (see Section S1.2, ESI †).Fig. 7a shows the results of multiscale coarse-graining for U CG (r) at different surface concentrations.The shown results imply that the effective interaction potential has a soft repulsive core, and an attractive tail begins to appear at smaller surface concentrations.Furthermore, at smaller surface concentrations, its dependency is weaker.A similar behavior for the effective interaction potential has been observed in our previous study 26 for surfactants with a different architecture adsorbed at a liquid-liquid interface.
To test the validity of the pairwise additivity assumption for CG interaction potentials and the accuracy of our force matching calculations, we have performed Monte Carlo simulations using the Metropolis scheme 40 for the corresponding two-dimensional coarse-grained systems.Fig. 7b shows a comparison between the radial distribution functions of (i) the projected surfactant centerof-masses obtained from the 3D reference simulation and (ii) the CG sites obtained from a 2D MC simulation, both for a system with a surface concentration of G = 0.118s À2 .The agreement is excellent, which validates our assumptions and force matching calculations.Excellent agreement between the pair correlation functions is also observed for the lower surface concentrations we examined (not shown).
The statistical mechanical foundation of the CG interaction potential defined by eqn (12) implies that U CG depends on the state point (T, G), as we confirmed by showing its dependency on G in Fig. 7a.One can perform 3D reference simulations for Fig. 6 A 3D snapshot, surfactant monolayer and surfactant center-of-masses projected onto the xy-plane for a system consisting of 600 HT 5 surfactants at a surface concentration of G = 0.118s À2 .Blue dots, yellow and red spheres stand for w, H and T particles respectively.Gray discs represent surfactant center-of-masses within the 2D interface, where a gas-liquid phase transition may occur.The right-hand-side graph shows the forcefield, its confidence intervals and the coarse-grained interaction potential (solid, dashed and dash-dotted lines respectively) obtained from multiscale coarse-graining.different surface concentrations and can construct U CG (r; G, T) over a range of concentrations and temperatures as we did in our previous study. 26Because the gas-liquid phase transition of Langmuir monolayers takes place at relatively low surface concentrations, and because we are exploring a new approach to determine the phase equilibrium of a model Langmuir monolayer for the first time, here we adopt the CG interaction potential evaluated at the lowest studied surface concentration and neglect its dependence on surface concentration and temperature in this relatively narrow parameter range.This way we can employ standard GEMC simulations and well-established liquid state theories.These approximations need to be carefully reviewed in future studies.

Gas-liquid phase transition of the monolayer
In this section, we investigate the 2D gas-liquid expanded phase transition of the monolayer for the CG interaction potential evaluated at G = 0.039/s 2 using an ''approximate'' density functional theory and ''exact'' Monte Carlo simulations in the Gibbs ensemble.

Gibbs ensemble Monte Carlo (GEMC)
The Gibbs ensemble technique 31 is a method to directly determine the phase coexistence properties of a system at fixed temperature from a single computer simulation with a relatively small number of E1000 particles.A system in the canonical ensemble is considered which consists of two periodic 2D simulation regions, one of them carrying the 2D liquid-expanded phase, the other the 2D gas phase, without a direct gas-liquid dividing line.These two subsystems can exchange particles and change their areas such that where N is the total number of particles and A is the total area, and the values with superscripts le and g correspond, respectively, to those in the liquid-expanded and gaseous subsystems.The equilibrium condition requires the temperature, chemical potential and surface pressure to be equal in the two coexisting phases, originally derived by J. W. Gibbs as 27,31 T le = T g , m le = m g and P le = P g .( 18) Three different trial moves are performed in the MC simulation of the combined system to satisfy coexistence conditions in eqn (18) subject to constraints eqn ( 17): (i) independent particle displacements within simulation regions le and g, (ii) area exchange and (iii) particle exchange between the two regions.While exhaustive details of the Monte Carlo simulations in the Gibbs ensemble can be found in ref. 27 and 31, we have collected the most important aspects together with technical details in the Section S2 (ESI †).Snapshots of the GEMC simulation of the 2D gaseous and liquid phases coexisting at T = 0.46e/k B and their corresponding pair correlation functions are shown in Fig. 8.The pair correlation function of the gas phase shows only one peak exactly at the same radial distance where the minimum happens for the CG interaction potential (Fig. 7a), whereas the liquid phase is more structured and one can identify second-nearest neighbor peaks.
The particle exchange move in the GEMC simulation provides us with the two chemical potentials m le and m g without any additional cost through where DU k ex is a test particle energy, L is the thermal de Broglie wavelength, and hÁ Á Ái k denotes the restricted expectation value for subsystem k.The surface pressures P le and P g in the two subsystem can also be obtained from the virial theorem as 41 We note that we ignore here a correction term to eqn (20) due to the density-dependence of U CG . 26The densities, surface pressures and chemical potentials of the two coexisting phases obtained from GEMC simulations will be compared with the results of DFT calculations in the following section.

Density-functional theory (DFT)
The coarse-grained sites representing the surfactant center-ofmasses interact within the 2D interface via a pairwise interaction potential U CG (r) which consists of both, short-range repulsive and attractive branches at large intermolecular distances.We use the Weeks-Chandler-Anderson scheme 42 and divide U CG into a repulsive reference part, U rep (r), and an attractive contribution, U att (r).They are shown in Fig. 9 and are defined by and where r min denotes the energetically preferred distance, where U min = U CG (r min ).We follow ref. 30 and define a family of intermediate interaction potentials via where increasing l from 0 to 1 corresponds to the gradual increase of the attractive tail and hence a linear change from a repulsive reference system characterized by U rep to the system of interest characterized by U CG .The 2D systems's intrinsic Helmholtz free energy F[r] is a functional of the system areal density r(r) and is exactly given by 30 with the Helmholtz free energy F 0 of the repulsive 2D reference system described by U 0 = U rep , while r 2 (r,r 0 ;l) is the pair areal density of the system in which particles interact with U l according to eqn (23).This pair density is not known in general for inhomogeneous systems, 43,44 but it is related to the radial distribution function g(r,r 0 ;l) through r 2 (r,r 0 ;l) = g(r,r 0 ;l)r(r)r(r 0 ) and within the random phase approximation (RPA) [43][44][45][46][47] given by r 2 (r,r 0 ;l) E r(r)r(r 0 ).
The RPA ignores all correlations between particles and is assumed to work quite well for weakly inhomogeneous systems such as gas-liquid interfaces [43][44][45][46][47] away from the critical region.
To calculate the free energy F 0 [r] of the repulsive 2D reference system, we use the Barker-Henderson (BK) scheme 48 which approximates F 0 [r] by thus employing the free energy areal density f hd (r) of a system of hard disks with temperature-dependent diameter d BH given by 48 d BH ðTÞ ¼ and visualized in Fig. S2 (ESI †).A local density approximation (LDA) is used in eqn (26) which is suitable for weakly inhomogeneous systems, e.g., for liquid-gas interfaces. 43,44We have calculated the effective hard-disc diameter at T = 0.9e/k B from eqn (27) as d BH = 1.816AE 0.009s which is of the same order but somewhat between the values obtained in Section 3.1 through fitting eqn ( 8) and ( 9) to the results of the P-A isotherm (Fig. 5).
To summarize, the Helmholtz free energy functional, F[r], after combining eqn ( 24)-( 26) can be approximated as 43,44 F½r % where all integrals extend over the domain of the 2D interface.
To calculate the required thermodynamic property f hd of the hard disc fluid we have used the hard disc equation of state 49 involving the packing fraction Z of hard discs.Using this equation of state, the chemical potential of hard disc fluid, m hd , is described as and the free energy density of the hard disc fluid, f hd = ÀP hd + rm hd , is obtained as This expression can be substituted into eqn (28) to estimate F 0 (r) of the repulsive 2D reference system.In the homogeneous limit of position-independent r(r), the Helmholtz free energy functional F[r] given by eqn ( 28) can be cast into the following form F[r] E Af(r) with the free energy areal density f (r) given by where a is determined by the available U att (r), The corresponding chemical potential and surface pressure in the uniform limit are given by mðrÞ ¼ df ðrÞ dr ¼ m hd ðrÞ À ar (34)   and Eqn (34) and (35) are the main relationships used next to determine the gas-liquid coexistence densities and the chemical potential.At a given temperature T the coexisting Fig. 9 CG interaction potential (dashed-line) and its attractive (dash dottedline) and repulsive branches (solid line) as obtained from the 3D reference simulation at G = 0.039/s 2 .
To estimate the critical point (r c ,T c ) from GEMC simulations, we fit the density difference of the coexisting phases to the scaling law 27,28,31,32 as where the critical temperature and the critical exponent are obtained as T c = 0.548 AE 0.005e/k B and b = 0.31 AE 0.01 respectively.Note that this exponent b is not equal to the corresponding critical exponent of the 2D Ising model and is in fact much closer to that of the 3D one.In order to estimate the critical exponent and the critical temperature more properly, knowledge of interaction potentials U CG (r; G, T) over the examined temperature and areal density ranges, and finite-size scaling analyses are necessary.We further fit the results of the GEMC simulations to the law of rectilinear diameters, 27,28,31,32 and estimate the value of the critical density as r c = 0.087 AE 0.001s À2 .The corresponding DFT estimates for the surface pressure and chemical potential at coexistence are plotted in Fig. 11 against the results of the GEMC simulations, obtained from ensemble averages described by eqn (19) and (20).The fact that the chemical potentials and pressures of the gaseous and liquid expanded phases in the GEMC simulations are equal within statistical errors validates the accuracy of our MC simulations and the sampling of equilibrium configurations.The DFT result (solid lines in Fig. 11b) for the chemical potential is in a good agreement with the result of the GEMC simulations.The agreement of the surface pressure (Fig. 11a) is very good at low temperatures T r 0.35e/k B but becomes increasingly inaccurate as T gets closer to the critical point.As is obvious from Fig. 10, and in part due to the approximations involved, the critical point resulting from the DFT calculation is located at a higher T and smaller r compared with the GEMC result.
Apart from calculating the areal densities, surface pressure, chemical potential, and free energy at coexistence, the form of the density profile across the gas-liquid expanded dividing line and the line tension can also be determined from DFT. [44][45][46][47] The grand potential of an inhomogeneous monolayer with density variations in the y-direction is given by where m = m(T) is the chemical potential at coexistence (obtained from solving eqn ( 36) and (37) simultaneously), and is constant throughout the system under the equilibrium condition.The equilibrium density profile r(r) fulfills the following variational condition 30,[43][44][45] dO½r drðrÞ ¼ 0: As a consequence, the density profile r(y) normal to the gasliquid expanded dividing line defining the x-axis must fulfill the following integral equation  where y and y 0 denote the transverse components of r and r 0 with respect to the 1D phase boundary.Eqn ( 42) is selfconsistently solved for r(y) by considering an initial density profile r (0) (y) and applying an iterative procedure [44][45][46][47] as until convergence has been reached, where r (i+1) and r (i) represent the i + 1 and ith iterations respectively (see Section S3.4,ESI †).The procedure is repeated for every T separately, and a profile r(y) obtained at a given temperature can serve as an initial r (0) (y) for a subsequent T.
The resulting density profiles across gas-liquid expanded dividing lines at different temperatures (Fig. 12) can be fitted very accurately to a hyperbolic tangent profile 51 where l 0 = l 0 (T) is the width of the dividing line region and y 0 represents the position of the equimolar dividing line.This procedure reveals that by increasing the temperature, the thickness increases and the dividing line region becomes more diffuse: while at T = 0.425e/k B the thickness is l 0 = 12.88 AE 0.14s, at the higher T = 0.502e/k B we obtain l 0 = 16.52 AE 0.09s within the DFT approach.
The grand potential of the inhomogeneous 2D system can now be determined by plugging the DFT equilibrium density profile, r(y), into eqn (40) and using eqn (28) as where, for example, Ð drf hd ðrðyÞÞ ¼ L Ð f fd ðrðyÞÞdy with L denoting the length of the dividing line, and N ¼ Ð drrðyÞ representing the total number of particles within the monolayer, which drops out from the DFT calculations.Now we have all ingredients to calculate the line tension of the gas liquid-expanded dividing line as 46 where A is the total area of the gas and liquid expanded phases together.The results of the gas-liquid line tension at different temperatures are shown in Fig. 13.The line tension decreases monotonically by increasing the temperature, 46 and should in principle vanish at the critical point where there is no distinction between liquid and gas phases. 50The results of the grand potential and line tension can also be used to investigate 2D gas liquid-expanded nucleation for the surfactant monolayer, see e.g.Zeng, 46 Zeng and Oxtoby. 47

Conclusions
Exploring G-LE phase transitions of Langmuir monolayers is extremely time-consuming via finite-size 3D molecular simulations.Starting out with MD simulations of a 3D reference model system containing a monolayer of insoluble surfactants at the vapor-liquid interface, we determined the effective lateral interaction potentials between surfactant center-of-masses at the interface.MC simulations of the resulting 2D systems show that these effective interaction potentials are quite successful in predicting pair correlation functions (and hence the structure) of surfactant center-of-masses projected onto the interface plane.At relatively small surface concentrations, the effective interactions exhibit an attractive tail suggesting that in principle a gas-liquid expanded phase transition can take place within the monolayer.For the effective 2D (simple) fluid of surfactant center-of-masses, described by these CG interaction potentials, we performed GEMC simulations and DFT calculations to investigate the 2D gas-liquid phase equilibrium, while these interactions could also be used to perform 2D simulations for large system sizes, and to explore large-scale self-assembly and domain formation within the monolayer.The DFT results for coexisting areal densities, surface pressure, and chemical potential are in semi-quantitative agreement with those of the GEMC simulations.
We have illustrated and applied these procedures for a system at a relatively small surface concentration, where we neglect the dependence of U CG on surface concentration G and temperature T. As a next step, it would be desirable to determine the CG interaction potential over a concentration range and construct  U CG (r; G, T) 26 at different temperatures.Then it is relatively straightforward to perform GEMC simulations for densitydependent interaction potentials, as discussed by Smith et al. 52 The complexity of the DFT calculations remains unchanged upon varying G and T. Such studies would then allow to evaluate the effect of the density-and temperature-dependence of U CG on the phase behavior, line tension, and density profiles.
Here we employed as a reference system a linear bead-spring model representing surfactants at the vapor-liquid interface of a monoatomic Lennard-Jones solvent. 8The presented scheme can be applied straightforwardly to other surfactant architectures such as branched phospholipids.Such effort could enable us to resolve the effect of molecular interactions and surfactant architecture on their phase behavior.

Fig. 1
Fig.1(a) An equilibrium snapshot of the 3D reference system with 600 HT 5 surfactants (300 per monolayer) at a surface concentration of G = 0.118s À2 .Blue dots represent w particles.Red and yellow spheres stand for H and T beads respectively.(b) Number densities as a function of z, normal to the 2D vapor-liquid interface; n w and n s denote w and surfactant center-of-mass number densities respectively, while n p stands for the total particle number density including w particles, H and T monomers.

Fig. 2
Fig. 2 Top view snapshots of HT 5 surfactants at different surface concentrations of (a) G = 0.0098, (b) 0.0196, (c) 0.0392 and (d) 0.0784s À2 ; remaining thermodynamic conditions as for Fig. 1.Yellow and red spheres stand for H and T beads, respectively.Grey dots represent w particles.

Fig. 4
Fig. 4 (a) Number density profiles of surfactants (based on the position of their center-of-masses) at different surface concentrations (symbols) and Gaussian fits (solid lines) for the 3D reference system.(b) Bond orientational order parameter for surfactants adsorbed at the monolayer as a function of surface coverage.

Fig. 7
Fig. 7 (a) Coarse-grained interaction potential at different surface concentrations G. (b) A comparison between pair correlation functions of (i) projected surfactant center-of-masses obtained from 3D reference simulations and (ii) the corresponding 2D coarse-grained system obtained from 2D Monte Carlo simulation, both at a surface concentration of G = 0.118/s 2 .

Fig. 8
Fig. 8 Pair correlation functions (a) of the 2D gas phase (blue line) and the liquid phase (red line) and the corresponding simulation snapshots of the liquid phase (b) with the density of r le = 0.1693s À2 and the gas phase (c) with the density of r g = 0.0190s À2 coexisting at T = 0.46e/k B .
3, ESI †) and then provide us with a unique m le (T) = m g (T) = m(r le (T)) etc. for comparison with the GEMC results for m le (T) and m g (T).

4. 3
Fig. 10 displays the phase diagram of the 2D fluid of surfactant center-of-masses obtained from DFT against the results of 2D GEMC simulations.The phase diagram consists of gaseous and liquid phases at temperatures above the critical temperature T c and their coexistence region below T c .It shows that DFT predicts semi-quantitative results for the coexisting areal densities with more accurate values for the vapor branch than the liquid branch.50To estimate the critical point (r c ,T c ) from GEMC simulations, we fit the density difference of the coexisting phases to the scaling law27,28,31,32 as

Fig. 10
Fig. 10 Phase diagram of the 2D fluid described by U CG in Fig. 9 obtained from DFT (solid line) and GEMC simulations (open circles).The left branches correspond to the gaseous phase and the right ones correspond to the liquid phase.The critical point (filled circle) is estimated by fitting GEMC results to the scaling law and the law of rectilinear diameters.

Fig. 11 (
Fig. 11 (a) GEMC results for the surface pressure in the liquid phase (open squares) and the gas phase (open circles) within the coexistence region below T c .The solid line shows the pressure at coexistence obtained from DFT for comparison.(b) GEMC results for the chemical potential in the liquid phase (open squares) and the gas phase (open circles).The solid line shows its value at coexistence obtained from DFT.

Fig. 12
Fig. 12 Density profiles across equilibrium gas-liquid dividing lines at different temperatures, obtained via DFT.The left and the right plateaus correspond to the 2D gas and liquid bulk densities, respectively.

Fig. 13
Fig. 13 The equilibrium gas-liquid line tension of the 2D fluid of surfactant center-of-masses, obtained via DFT, as a function of temperature.The line tension vanishes at the critical T c predicted by DFT.