Application of molecular simulations: Insight into liquid bridging and jetting phenomena

Molecular dynamics simulations have been performed on pure liquid water, aqueous solutions of sodium chloride, and polymer solutions exposed to a strong external electric field with the goal to gain molecular insight into the structural response to the field. Several simulation methodologies have been used to elucidate the molecular mechanisms of the processes leading to the formation of liquid bridges and jets (in the production of nanofibers). It is shown that in the established nanoscale structures, the molecules form a chain with their dipole moments oriented parallel to the applied field throughout the entire sample volume. The presence of ions may disturb this structure leading to its ultimate disintegration into droplets; the concentration dependence of the threshold field required to stabilize a liquid column has been determined. Conformational changes of the polymer in the jetting process have also been observed.


Introduction
With advance in computer technology, the development of efficient simulation techniques, and the availability of commercial or freeware packages, molecular simulation methods [1,2] have become the main, and in some instances the only tool to study fluid systems at the molecular level. A primary output of molecular simulations (and target of theories) have typically been thermophysical properties (equation of state, heat capacity, viscosity, etc.), entering macroscopic theories of the studied systems, and their relation to various kinds of underlying intermolecular interactions. However, they make it also possible to study various processes and this approach has become thus commonplace in many other fields, including mechanical engineering and material design, biology, and even medicine, often shedding light on problems that are inaccessible by other methods. Examples include the rheology of lubricants [3,4] and protein self-assembly [5], potentially related to Alzheimer's and Parkinson's diseases. In these applications, the simulations need not necessarily yield numerical values of physical quantities, but they may help to qualitatively understand and explain the underlying molecular mechanism of the observed phenomena.
A similar type of problems whose molecular nature is not yet fully understood involves liquid surfaces exposed to strong external electric fields, their disturbance and subsequent development. These problems are encountered in a variety of fields, e.g., physics of the atmosphere (ice formation in supercooled water), biology (transport through biological cell membranes), and mineralogy (reactivity at the mineral-water interface). Interesting examples of processes in which the external field is the driving force are the floating water bridge (FWB), whose application potential has not been explored yet [6], and the production of nanofibers using needleless electrospinning technology [7]. The macroscopic floating bridge arises between two beakers filled with a liquid (not necessarily water) and with immersed electrodes [8], see figure 1 (A). The beakers are initially in contact and a free hanging structure is formed when they are then slowly pulled apart. A nanoscale analogue to the floating bridge (a nanobridge) is observed between the tip and sample surface in atomic force microscopy (AFM) [9]. Several scattering experiments using different techniques have been used to elucidate the microscopic mechanism which holds the FWB together but their results were inconclusive or even contradicting [10][11][12][13]. And, as for molecular simulation studies, we are aware of only two attempts. Skinner et al. [13] used the simulation setup not fully representing the FWB but only its bulk-like interior and did not find any preferential orientation of the water molecules in the bridge. Cramer et al. [14] studied a drop attached to a surface and observed its deformation into a column-like structure, though never fully developing a bridge. Their main result was the finding that whatever strength of the field was used, the hydrogen bond network of the water molecules remained practically intact.
In needleless electrospinning, the external field destabilizes the surface of the liquid. It gives rise to the formation of a Taylor cone which ultimately looses stability and an eruption of a jet is observed; these stages are schematically shown in figure 2. The photograph in figure 1 (B) depicts the experimental observation of jetting. After flying through a short linear section, the jet undergoes the so-called whipping stage and the resulting fiber is then collected at the other electrode. With the exception of a recent paper [15], we are not aware of any other atomistic simulation study dealing directly with the above described process. Other available simulations related to this process focus on the structure formation of polymer chains [16], the effect of solvent on the produced nanofibers [17], and on spraying ionic liquids [18].
In a recent paper [15] we made the first attempt to apply molecular simulations to the electrospinning process to examine (i) to what extent simulation can be applied to such a dynamic process and (ii) what results it can provide. We used three different methodologies, all of which yielded qualitatively the same results. We have continued in this research in which, in addition to the qualitative pictures, the gain of also some quantitative characterizations has been attempted and the obtained results are reported in this paper. We consider both the solution of sodium chloride and a model polymer solution, specifically the solution of polyethylene glycol (PEG). For completeness and comparison, also pure water has been included in this study. In the next section we provide the main technical details of the simulations, and in section 3 we present the obtained results. These include the stability analysis in the dependence on the strength of the field, an analysis of the orientational arrangement of the water molecule dipoles and the effect of the presence of ions thereupon in dependence on the electrolyte concentration, changes in the polymer conformation during the process, and relevant kinetic properties.

Methodology
For molecular dynamics simulations reported in this paper, we used the GROMACS package (versions 4.5 and 4.6) [19]. A leap-frog integrator with time steps ranging from 1 to 2 fs was employed and the particle-mesh Ewald method [20] with cubic interpolation and tinfoil boundary conditions were used to treat electrostatic forces whenever periodic boundary conditions were applied. The technical parameters of the Ewald summation (cutoff, r c , separation parameter, α, grid spacing, etc.) were being adjusted for different setups individually to optimize computational cost while maintaining reasonable accuracy. The relative strength of the direct potential at the cutoff (GROMACS ewald-rtol parameter), i.e., erfc(αr c ), was always set to the value 10 −5 .
To keep contact with the previous study [15], for simulations with pure water and sodium chloride solutions, we used the TIP3P water model [21] and the associated Joung and Cheatham [22] ionic force fields. The solution of polyethylene glycol (PEG) was modelled by the GROMOS 53A6 OXY+D [23] force field with SPC water.
In order to approach the proper scale of the investigated processes as closely as possible, and to make simultaneously the simulation feasible using the available computing facilities, the number of molecules in the simulations ranged from 3400 up to 55000. In simulations with polymer, two 10.5 kDa PEG chains of the formula HO-[-CH 2 -CH 2 -O-] n -H, n = 239, were placed into the simulation box along with water and ions.
The supporting square underlay in the simulations without periodic boundary conditions was realized by a square grid of Lennard-Jones (LJ) atoms with fixed positions. The spacing between neighbours in the grid was 0.14142 nm. The contact angle depends on the mutual interaction between the droplet molecules and the underlay. We fixed, therefore, the parameters to the values corresponding to a weakly hydrophobic material to provide a sufficiently firm attachment while preventing the undesired spreading of liquid over the underlay surface. Furthermore, in order to prevent the droplet from moving over the edge of the underlay, a repulsive rim made of sites with the ∼ r −12 repulsive potential was added. The parameters of the LJ potential C 12 r −12 − C 6 r −6 were, respectively, C 12 = 1.234 × 10 −6 kJ mol −1 nm 12 and C 6 = 5.852 × 10 −4 kJ mol −1 nm 6 for the underlay, and C 12 = 4.937 × 10 −6 kJ mol −1 nm 12 for the rim.
For the interaction of the underlay with the liquid, geometric mean combining rules were applied to both constants (GROMOS default).
The intensity of the applied electric field varied within the range from 0.01 to 2 V nm −1 . The use of a field of such strength is consistent with simulation studies on water reported in the literature [14,24,25]. In real electrospinning experiments, the applied voltage is, typically, between 10 and 30 kV. The exact strength around the apex of the Taylor cone depends on the electrode separation and the geometry of the apex. A rough estimate gives values between 1 and 100 V nm −1 which also provides support for the strength of the field used in the simulations.
Before the field was applied, all systems were preequilibrated in an NV T ensemble at 298.15 K using the stochastic velocity rescaling method of Bussi et al. [26] with correlation time 0.1 ps. In simulations involving jet formation, no thermostat was in use after the field had been turned on, while in bridge stability studies, a thermostat was applied during the entire simulation. Ions were excluded from velocity scaling to avoid interference with their drift.

Results and discussion
A typical development of the jet observed in simulations with dilute solutions of NaCl is shown in figure 3. Qualitatively the same behavior was also observed for pure water. Without the external field, the solution remains on the underlay in the form of a droplet, figure 3 (A). When the homogeneous external field in the direction perpendicular to the underlay is switched on, the droplet gradually develops into a cone, figure 3 (B), and then it forms a pillar resembling a jet in electrospinning. Apart from a small number of individual water molecules belonging to its vapor phase, the jet remained continuous. It means that no molecular clusters broken off the jet (which is typical of electrospraying) were observed. Since the results obtained for pure water and a dilute electrolyte were qualitatively the same, one can deduce that ions, i.e., charge bearing particles, are not the driving force behind the origin of jetting. One can thus speculate that the building blocks of the jet are the dipoles (either permanent or induced) of molecules. That the dipoles do form a chain aligned along the external field has been confirmed and this is demonstrated below. Nonetheless, the internal structure of pure water may be more complex because even at a strong external field the hydrogen bonding network remains intact despite the alignment of a dipole as shown by Cramer et al. [14].
In figure 4 we show the same sequence as above in the case of a more concentrated solution. At the beginning, the droplet follows the same evolution pattern as discussed above but at a later instant, instead of the development into a stable jet, small clusters formed by a central ion (hydrated ion) leave the bulk structure and accelerate away along the field lines of the external field. Although reality may be more complex, it is evident that the ions do disturb the structure observed in the case of pure fluids. The interaction between the ions and water molecules (hydration) may overcome the effects of external field on water (polarization). This intuitive deduction is confirmed in figure 5. The strong polarization of water molecules and their alignment along the field in the jet is demonstrated in panel (A). Panel (B) reveals that the ion Na + makes the water molecules in its vicinity reorient with their oxygen atom to  be closer to it. Similarly, the Cl − ion has the same effect with the difference that now it is the hydrogen atom which is pulled closer to the ion (C). This reorientation of water molecules disturbs thus the original alignment of the water dipoles which in the case of sufficiently high concentration leads to disintegration of the column-like structure and spraying or collapse of the bridge.
In our feasibility study [15] we considered only pure water and electrolyte solutions. In real electrospinning, it is polymer which in the end forms the nanofiber. Another step towards real electrospinning must also involve polymer molecules; we have also performed, therefore, simulations on a polymer solution. It is known that the majority of simulations with polymer solutions are carried out either on lattice or considering the solvent only implicitly. This is not suitable for studying the molecular nature of the electrospinning process and there are two questions to be answered: (i) whether a picture of the pro- cess can be obtained within a reasonable CPU time using explicit solvent and (ii) how large the polymer molecule should be to make simulations feasible and yet realistic. In figure 6 we show jet development from the polymer solution which does not differ much from the scenarios discussed above. The polymer molecule changes its conformation along with deformation of the droplet, figure 6 (B). The polymer chain uncoils when it is forced to enter the jet, figure 6 (C). A detailed investigation of the changes of the polymer conformations in dependence on the polymer composition, its structure, the strength of the field, electrolyte concentration, etc., is beyond the scope of this study and deserves an independent investigation. Nonetheless, the presented results show that such a study with polymer solution is nowadays feasible.
In figure 7 we show the results on the dynamics of the onset of jetting: the average mass flow (MF) in its dependence on the distance from the axis running through the jet's center-of-mass and parallel to the external field and also on the distance from the underlay. As it is seen, (A), deep in the bulk solution and far from the bottom of the jet, the effect of the external field on MF z is negligible. However, as soon as the solution enters the region of the jet, it sharply accelerates in the direction of the field. In figure 7 (B) we show the motion of the solution in the direction towards the jet's axis. Despite a large statistical noise, we see that on the surface of the solution in the vicinity of the base of the jet there is a region where the solution is more strongly drawn into the jet. We may thus gather that a relatively large part of the jet mass originates from the surface layer around the jet.
All results commented on above correspond to situations where the external electric field is strong enough to initiate transformation of the liquid surface into a jet. If the field was suddenly weakened under the jetting threshold at any instant during the simulation, the jet would eventually vanish. However, one can use quite a different setup, where we anchor a liquid pillar on both ends, forming a liquid bridge. In contrast to a jet, being a dynamical phenomenon, the bridge can exhibit a thermodynamic (meta)stability (in the absence of drifting ions) and can exist under much a weaker field than that needed for its build-up from the liquid surface. An illustration of this phenomenon is given in figure 8 showing snapshots from the simulation of the NaCl aqueous solution, which initially forms a column continuous in the z-direction, i.e., in the direction of the applied field, through periodic boundary conditions, see panel (A). This setup elegantly avoids attaching the ends of the bridge to a particular material by simulating only the middle portion of the bridge self-anchored by its periodic images. The bridge was originally formed from a liquid slab by a strong field perpendicular to the liquid surface, in a way similar to all the above-described simulations. Subsequently, a column was subjected to a series of weakened fields, while the stability of the structure, i.e., whether a column remains continuous or breaks, was assessed. When the field was weak enough, increasingly prominent bulging appeared along the column, see figure 8 (B), resulting in  Figure 9. Dependence of the electric field strength necessary for a stable bridge of aqueous NaCl solution to be formed on its concentration. All the initial bridges were 4 nm in diameter. Circles denote the stable bridge and crosses mark its breakup. the ultimate disintegration of the continuous liquid column into droplets, see panel (C). Figure 9 maps the stability of the established bridge with respect to the strength of the field and the electrolyte concentration. In accord with the above-given explanation of the effect of ions on the jet stability (discussion pertaining to figure 4), one can see that as concentration increases, the critical intensity of the electric field needed to maintain a stable bridge also increases. It is to be expected that this behavior also depends on the dimensions of the bridge; however, this was not investigated within the present study. All the bridges used in the stability assessment had diameters around 4 nm. The resulting field intensities are substantially larger than those usually applied in experiments with the FWB of macroscopic scales. One can thus expect that the critical field strength will decrease with increasing bridge diameter. However, a proper molecular study of this trend will require further extensive simulations.

Conclusions
Our study has shown that it is possible to model complex field-induced processes at liquid surfaces, such as the formation of the electrospinning jet or the collapse of the water bridge, at the atomistic level by molecular simulation. We have shown how such simulations can be performed and have given several examples of the microscopic phenomena which are experimentally inaccessible whilst being important for proper understanding of the behavior of aqueous solutions in electric fields. In contrast to intuitive views on the mechanism of the origin of jetting and bridging phenomena, our simulations show that the essential building blocks of the liquid column are the water molecules, and that this is qualitatively independent of the presence of ions in solutions. The ions, such as Na + and Cl − , destabilize both the growing jet and the water bridge. The primary mechanism leading to the destabilization is the disruption of the water orientational polarization in the ions' vicinity. The stability thus decreases with growing ionic concentration. As regards the dynamics of jetting, the simulations indicate that the particles are drawn into the jet preferentially from the liquid surface near the base of the jet and subsequently accelerate in the direction of the field. Each of the above conclusions, however, provokes new questions, which can be reliably answered only after further extensive simulations are performed as well as thorough analysis.