A Molecular-Level Picture of Electrospinning

Abstract

Electrospinning is a modern and versatile method of producing nanofibers from polymer solutions or melts by the action of strong electric fields. The complex, multiscale nature of the process hinders its theoretical understanding, especially at the molecular level. The present article aims to contribute to the fundamental picture of the process by the molecular modeling of its nanoscale analogue and complements the picture by laboratory experiments at macroscale. Special attention is given to how the process is influenced by ions. Molecular dynamics (MD) is employed to model the time evolution of a nanodroplet of aqueous poly(ethylene glycol) (PEG) solution on a solid surface in a strong electric field. Two molecular weights of PEG are used, each in 12 aqueous solutions differing by the weight fraction of the polymer and the concentration of added NaCl. Various structural and dynamic quantities are monitored in production trajectories to characterize important features of the process and the effect of ions on it. Complementary experiments are carried out with macroscopic droplets of compositions similar to those used in MD. The behavior of droplets in a strong electric field is monitored using an oscilloscopic method and high-speed camera recording. Oscilloscopic records of voltage and current are used to determine the characteristic onset times of the instability of the meniscus as the times of the first discharge. The results of simulations indicate that, at the molecular level, the process is primarily driven by polarization forces and the role of ionic charge is only minor. Ions enhance the evaporation of solvent and the transport of polymer into the jet. Experimentally measured instability onset times weakly decrease with increasing ionic concentration in solutions with low polymer content. High-speed photography coupled with oscilloscopic measurement shows that the measured instability onset corresponds to the formation of a sharp tip of the Taylor cone. Molecular-scale and macroscopic views of the process are confronted, and challenges for their reconciliation are presented as a route to a true understanding of electrospinning.

1. Introduction

Electrospinning is a technique utilizing high-voltage (HV) electric fields to produce nanofibrous materials from polymer solutions or, less commonly, polymer melts [1]. In the electrospinning device (for a schematics see Figure 1), the strong electric field destabilizes the surface of the liquid connected to a terminal of a high-voltage power supply, and the liquid subsequently forms a pointed meniscus, the Taylor cone. A discharge in the form of a thin liquid jet then emanates from the tip of the cone. The liquid jet is not always stable and may disintegrate into liquid aerosol leading to an electrospray, which is itself an important phenomenon in applied physics [2,3]. When the jet arrives at the counter-electrode, a collector, it deposits as a solid nanofiber. In the case of solution spinning (which is the most common one), solidification is facilitated by rapid evaporation of the solvent. Bending, or whipping, instability occurring after a short linear path of the jet significantly contributes to elongation and thinning of the jet [4,5]. By depositing many layers of nanofibers, a non-woven fabric is created on the collector showing exceptional porosity and surface-area-to-volume ratio [6]. Owing to their unique properties and versatility of their preparation, electrospun nanofabrics find a wide variety of applications, such as air filtration [7,8], oil/water separation [9], drug delivery [10,11], or tissue engineering [12]. With the outbreak of the novel coronavirus pandemic (COVID-19) in 2020, electrospun materials have gained attention as a promising active layer for protective face masks [13].

The theoretical description of electrospinning and related phenomena has predominantly relied on electrohydrodynamic models formulated within continuum mechanics [5,15,16]. The fundamentals of theoretical understanding of the field-induced phenomena at liquid interfaces related to electrospray and electrospinning were laid by Lord Rayleigh [17], John Zeleny [18], and most notably sir Geoffrey Ingram Taylor [19,20,21], who studied a conical formation occurring on liquid surfaces in strong electric fields, the formation which today bears his name: the Taylor cone. Conventional electrohydrodynamic approach to weakly conducting fluids in electric fields, the Taylor–Melcher leaky dielectric model [16], incorporates both dielectric polarization and ionic conduction in terms of the Maxwell stress tensor and conservation equations. If the charge relaxation time, ${\tau }_{\mathrm{c}}$, is much lower than the hydrodynamic relaxation time, ${\tau }_{\mathrm{h}}$, the behavior of a leaky dielectric liquid approaches that of a perfect conductor, and when the opposite is true, i.e., ${\tau }_{\mathrm{c}}>{\tau }_{\mathrm{h}}$, the liquid behaves as an insulator. The former situation is assumed to take place when the Taylor cone builds, whereas the latter is believed to exist in a flying jet [5]. For a modern review on continuum modeling of electrospinning, we refer the reader to the monograph of Ko and Wan [22] or the chapter by Rafiei [23].

Apart from the conventional continuum modeling one can also find reports on methods involving discrete macro- or mesoscopic elements. The entire jet can be modeled as a chain of charged particles connected by visco-elastic links with one another as, e.g., in works of Reneker et al. [4] or Lauricella et al. [24,25]. Mesoscopic modeling incorporating individual polymer chains was performed by Liu et al. [26,27,28,29], who employed the dissipative particle dynamics to model polymer melt electrospinning using a coarse-grained model of chains. The dissipative particle dynamics was also applied by Joulaian et al. [30] to model the initial stages of the formation of the Taylor cone.

Truly atomistic simulations of electrospinning with explicit solvent are scarce in the literature. Jirsák et al. published several works [31,32,33,34,35] developing and applying an approach to model electrospinning by atomistic molecular simulations in several geometries of the simulation cell using the rigid water model, simple sodium chloride ions and model poly(ethylene glycol) (PEG) chains. Jirsák et al. concluded that even a perfect dielectric (model water) was able to jet in electric field and revealed a link between liquid jetting and bridging. In a similar, independent study, Wang et al. [36] investigated the behavior of a polymer solution droplet in the electric field and the influence of ions on spinnability.

The above-mentioned previous works [31,32,33,34,35] co-authored by the first author represent a continuing effort to tackle the atomistic modeling of electrospinning methodologically and to understand the process on a fundamental, molecular level. In the first article [31], Jirsák et al. partitioned the process of electrospinning into zones accessible to modeling (free liquid surface, the tip of the Taylor cone, and the linear segment of the liquid jet) and proposed several strategies for implementation of molecular simulation in the respective zones. Application of the developed methodology to pure water and aqueous sodium chloride revealed that at atomic scale ions play only a marginal role in jet initiation; nevertheless, their presence destabilizes the jet. In the subsequent paper [32], the authors noticed an interesting analogy between the liquid jet in electrospinning and a cylindrical structure created by water in electric field, the ‘floating water bridge’ [37]. In paper [33], the authors further investigated the influence of ions on the stability of field-induced liquid cylinders, monitored the mass flow in jets and added a model of polymer to the system (PEG). In Ref. [34], Jirsák et al. analyzed polymer conformation and solvation during jetting. They concluded that polymer chains unfolded when entering the jet and that they remained well solvated despite progressing fragmentation of the solvent. The invited chapter [35] presented an comprehensive overview of the developed methodology and summarized the achieved results.

In the present contribution, the authors work further towards understanding the molecular mechanism of electrospinning, especially its initial stages: formation of the Taylor cone and the onset of jetting. Following the methodology developed earlier [31,32,33,34,35], systematic simulations of saline aqueous solutions of PEG are performed and analyzed with respect to a wide range of dynamic and structural properties. Unlike the previous studies, the present work adds laboratory experiments performed on the solutions of composition comparable to that used in simulations. An oscilloscopic method developed by P. Pokorný and coworkers at the Technical University of Liberec [38,39] is used to measure the latency time before the first Taylor cone builds and discharges.

The remaining part of the article is organized as follows. In the next section, methods utilized in molecular simulations and experiments are described along with all parameters necessary for reproduction. The following section, Results and Discussion, compiles the data produced by simulations and experiments, offers the interpretation of the data, and finally discusses how molecular and macroscopic observations might be related. Limitations of the present study are also discussed. The last section, Conclusions, summarizes main findings of the work and outlines challenges for future research.

2. Methods

2.1. Computational

2.1.1. Molecular Simulations

Classical molecular dynamics simulations of a nanodroplet of polymer solution on a solid underlay in strong electric fields were performed. GROMACS 4 package [40,41] was employed for all simulations, using the GROMOS 53A6${}_{\mathrm{OXY}+\mathrm{D}}$ force field [42] with simple point-charge (SPC) water model [43]. The leap-frog integrator with time step 2 fs was used. Bond lengths in polymer chains were constrained by P-LINCS [44] algorithm, SETTLE [45] was used to constrain a rigid water geometry.

The simulated polymer solutions comprised of three components: water, sodium chloride, and hydroxyl-terminated poly(ethylene glycol). Two different numbers of –O–CH${}_2$–CH${}_2$– monomer units, 119 and 239, were used, corresponding to PEG molecular weights 5260 Da and 10,547 Da, respectively. The number of water molecules in all solutions was around 23,700, whereas the population of other components varied. For sodium chloride, $N_{\mathrm{NaCl}}=0$, 18, 36, or 72 of Na${}^{+}$ Cl${}^{-}$ ion pairs were used. For 5260-Da polymer, $N_{\mathrm{PEG}}=8$, 12, or 16 chains were used, while for 10,547-Da polymer, $N_{\mathrm{PEG}}=4$, 6, or 8 chains were used. Initial configurations were prepared in a cubic cell using the highest concentrations of ions and polymer. The cube of the solution was then left to adhere to the surface of an underlay during a pre-equilibration stage until it formed a hemispherical droplet. Samples with lower concentrations were built by deleting excess solute molecules from the initial droplet of concentrated solution. The initial configurations resulting from the pre-equilibration step were then used for equilibration described later.

The support of a droplet was formed by a planar disc, ca. 20 nm in diameter, made of 15,384 Lennard-Jones (LJ) sites with positions fixed onto a square grid at $z=0$ with spacing 0.14142 nm. The LJ parameters of the underlay sites, $C6=5.851561\times {10}^{-4}$ kJ mol${}^{-1}$ nm${}^6$ and $C12=1.234321\times {10}^{-6}$ kJ mol${}^{-1}$ nm${}^{12}$, were adjusted to maintain a reasonable contact angle of the droplet close to 90${}^{\circ }$. Geometric combination rules were applied for cross-interactions of the underlay with other components. In order to prevent accidental ‘spilling’ of the liquid over the rim to the other side of the disc, the parameters of 912 sites of the rim (two outermost circles) of the underlay were modified to render a purely repulsive interaction by zeroing $C6$, keeping $C12$ unchanged. See Figure 2 for a graphical representation of the underlay. Figure 3 shows the underlay with a droplet.

A typical simulation protocol involved (i) a long 10-ns equilibration under periodic boundary conditions (PBC) in all directions at constant temperature $T=298.15$ K maintained by the Berendsen thermostat [47] with time constant $\tau =0.25$ ps, (ii) a short 1-ns equilibration without PBC using the stochastic velocity-rescaling thermostat of Bussi et al. [48] with the same T and $\tau$, and finally (iii) a non-thermostated production run without PBC under a stationary and uniform electric field applied in the normal direction to the underlay plane (positive z-axis). The production run lasted tens to hundreds of picoseconds. Between stages (i) and (ii), water molecules with oxygens located farther than 10 nm from the center of the disc were deleted to limit the number of molecules initially present in vapor phase.

In simulations with PBC, a cubic simulation box with edge length 25 nm was used. Electrostatic forces in an infinite periodic system were treated by the particle-mesh Ewald method [49] with cut-off 2 nm, grid spacing 0.3 nm, cubic interpolation, and the relative strength of the direct potential ${10}^{-5}$ at cut-off [41]. Lennard-Jones forces were neglected beyond a cut-off distance of 1.5 nm. In production runs without PBC, all site–site interactions were included in the calculation of the force, i.e., there was no potential cut-off applied.

Most of the production simulations were done with the external field intensity of 1.5 V/nm, but the effect of the strength of the field was also studied, using a series of simulation under 1.0, 1.5, 2.0, and 2.5 V/nm for selected systems. For each composition and field strength used, a single production trajectory was generated.

The national network of computer clusters Metacentrum [50] was used for all simulations, and typical computation on 32 CPU cores (AMD Opteron 6274 at 2.5 GHz) took over 11 days to equilibrate the system and another 14 hours to generate 100 ps of production trajectory.

2.1.2. Analysis of Trajectories

During the production runs, configurations of all sites were stored after each 1000 integration steps, i.e., every 2 ps, for subsequent analysis. Precision of the saved configurations was ${10}^{-4}$ nm. For selected systems also velocities were stored to determine the kinetic energy of ions.

The cumulative net mass transfer through a plane at $z=20$ nm was calculated as the total mass of atoms with $z>20$ nm by the given time instant minus the mass of atoms with $z>20$ nm at the start of production (to correct for the offset caused by initially present water vapor). The time of the jetting onset was determined as the time when the cumulative mass transfer through $z=20$ nm overcame the threshold of 2000 Da. The mass transfer through $z=20$ nm of PEG chains alone was also calculated and for the corresponding onset time the threshold of 500 Da was used.

The z-position of the center of mass of all PEG chains, $z_{\mathrm{PEG}}$, was numerically differentiated with respect to time to estimate the z-component of the velocity, $v_z$, of the center of mass of PEG at time instant t, i.e.,

$$v_z\left(t\right)\approx \frac{z_{\mathrm{PEG}}(t+\Delta t)-z_{\mathrm{PEG}}\left(t\right)}{\Delta t},$$

(1)

where $\Delta t=2$ ps is the trajectory sampling timestep.

The evolution of the dipole moment of the whole system was calculated during production run. Two components of the dipole moment were monitored: (i) the ‘dielectric’ dipole moment of water and PEG and (ii) the ‘ionic’ dipole moment arising from sodium cations and chloride anions.

The evaporation rate of water solvent was estimated based on counting one-particle water clusters at different time instants along the trajectory. Simple criterion of oxygen–oxygen distance less than 0.35 nm was used to define a cluster bond between particles. Water molecules not having any other water molecules bound to them by a cluster bond were considered to be a part of the vapor phase. The average evaporation rate in a given time interval was determined as the change of the number of vapor molecules divided by the length of the time interval. For all samples, three evaporation rates were systematically monitored: the average rate in the first 36 ps of production, the average rate in the second 36 ps of production, and the average rate in the first 72 ps of production (i.e., a mean value of the two preceding 36-ps rates).

In order to characterize the effect of ions on the evaporation of solvent, an attempt was made to evaluate the amount of energy dissipated from the ions (accelerated by the field) to the rest of the system. The dissipated ionic energy, $E_{\mathrm{dis}}$, was defined as the loss of the energy of ions between two chosen time instants $t_1<t_2$:

$$E_{\mathrm{dis}}(t_1,t_2)=E_{\mathrm{ion}}\left(t_1\right)-E_{\mathrm{ion}}\left(t_2\right).$$

(2)

The total ionic energy $E_{\mathrm{ion}}$ includes the kinetic energy of ions and the potential energy of both ion–field and ion–ion interactions:

$$E_{\mathrm{ion}}\left(t\right)=\frac{1}{2}\sum _i^{\mathrm{ions}}{m}_i{v}_i^2+\sum _i^{\mathrm{ions}}{q}_i\phi \left({\overrightarrow{r}}_i\right)+\frac{1}{4\pi {\epsilon }_0}\sum _{i<j}^{\mathrm{ions}}\frac{q_i{q}_j}{|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j|},$$

(3)

where $m_i$ is the mass of the i-th ion, $v_i$ its velocity, $q_i$ its charge, and ${\overrightarrow{r}}_i$ its position vector. For electrostatic potential $\phi$ it holds $\overrightarrow{E}=-\nabla \phi$. Therefore, with $\overrightarrow{E}=(0,0,E_z)$ we have

$$\phi (x,y,z)=-z{E}_z+\mathrm{const}.,$$

(4)

where the term arising from the constant vanishes in Equation (3) due to electroneutrality.

2.2. Experimental

2.2.1. Materials and Solution Preparation

The following chemicals were used for the preparation of solutions:

• Ultrapure water (18.2 M$\mathrm{\Omega }$ cm) was prepared by Purelab Flex 2 system (Veolia Water Systems, Ltd., High Wycombe, UK).
• Poly(ethylene glycol), ${\overline{M}}_{\mathrm{n}}$ = 6000 Da, was purchased from Sigma-Aldrich (Prague, Czech Republic).
• Sodium chloride p. a. was purchased from Penta (Prague, Czech Republic).

A total of 12 solutions were prepared using several different concentrations of PEG and NaCl in water. All solutions were based on ca. 100 g of water. The 12 solutions resulted from the combination of three different amounts of added PEG: 9.9, 14.8, and 19.7 g (denoted as a, b, and c PEG series, respectively), and four different amounts of added NaCl: 0.00, 0.25, 0.49, and 0.99 g (0, 1, 2, and 3 NaCl series). The above masses are only approximate, actual masses were determined to four or five significant digits for each solution. For the precise concentrations of PEG and NaCl in the prepared solutions, see Results and Discussion. Conductivities of all prepared solutions were measured by a conductivity meter (COND 51+ with 2301T cell, XS Instruments, Carpi, MO, Italy).

2.2.2. Instrumentation

In order to monitor the instability onset, a circuit was built as shown on Figure 4. A steel working electrode was connected, through a fast-acting mechanical switch, to the positive terminal of a high-voltage direct-current power supply (AU-60P0.5-L, Matsusada Precision, Inc., Kusatsu, Japan), while the ground terminal of the supply was connected to a plate-like collector electrode. A digital oscilloscope (DS1102CD, RIGOL Technologies, Inc., Munich, Germany) was used to monitor voltage and current through the electrodes. Working electrode with a droplet was recorded by a high-speed camera (i-SPEED 720, Olympus Europa SE & Co. KG, Hamburg, Germany) with framerate 20,000 fps. For high-speed camera recording, the working electrode was illuminated by a light source (ILP-1, Olympus Europa SE & Co. KG, Hamburg, Germany) with color temperature 8200 K. The wiring of the oscilloscope to the circuit, see Figure 4, allowed for monitoring electrode voltage on channel 1 (through R1–R2 voltage divider) and electrode current directly proportional to voltage on channel 2. The delay of the rise of the current signal behind the rise of the voltage signal (at closing the switch) was considered the time of the first discharge marking the onset of the instability of the meniscus.

Two slightly different experimental setups I and II, were used, see

Table 1. Technical parameters of experimental measurements.

Parameter	Experimental Setup
	I	II
working electrode diameter	10 mm	1.1 mm
working electrode axis	vertical	horizontal
est. droplet volume	0.2 mL	0.5 $\mathsf{\mu }$L
collector shape	disc	square
collector dimensions	⌀ 120 mm	$300\times 300$ mm${}^2$
collector material	aluminium	steel
electrode separation	70 mm	100 mm
voltage	20 kV	10 kV
current limit	0.15 mA	0.10 mA

. Setup I was used for both voltage–current delay measurement and high-speed camera recording, whereas setup II was used only for the high-speed camera recording. Figure 5 shows droplets placed on working electrodes in both setups. The high-speed camera was coupled with multifunction IO device (USB-6216-BNC, National Instruments, Austin, TX, USA) with two analog inputs (AI0 and AI1) connected to the same points as the channels of the oscilloscope in order to acquire voltage and current data synchronized with camera recording (5 data samples per frame).

Experiments were conducted according to the following protocol:

• The working electrode was briefly connected to the ground terminal to neutralize any residual charge.
• A droplet of polymer solution was placed onto the working electrode by a pipette.
• With switch S1 open, the HV power supply connected to the circuit was turned on.
• Switch S1 was closed, voltage and current signals were recorded by the oscilloscope. The behavior of the droplet was captured by the high-speed camera.
• HV supply was turned off. The working electrode was again discharged with a grounded wire.

In the case of voltage–current delay measurement, all experiments were repeated a minimum of 10 times. The droplet was placed so that it covered the whole upper surface of the rod electrode. The experiments where the droplet spilled over the edge before or during measurement were discarded. After each measurement, the electrode was carefully cleaned of the solution before placing a new droplet. Solutions with lower concentrations of ions were always placed prior to the more concentrated ones to minimize the influence of ions left adsorbed on the electrode surface.

2.2.3. Analysis of Data

Oscilloscopic records were analyzed directly in the oscilloscope. The time difference between the crests of the rising edges of voltage (channel 1) and current (channel 2) signals was considered to be the time of the instability onset, see Figure 6. High-speed camera recordings were processed by i-SPEED Software Suite (Olympus) [52].

3. Results and Discussion

3.1. Computational

3.1.1. Composition of Systems

Compositions of all systems used for MD simulations are summarized in

Table 2. Numbers of water molecules in production runs of all the simulated systems.

${\mathit{N}}_{\mathbf{PEG}}$	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36	72
$n=119,M_{\mathrm{r}}=5260$
8	23,684	23,682	23,684	23,695
12	23,687	23,695	23,693	23,707
16	23,692	23,682	23,686	23,703
$n=239,M_{\mathrm{r}}=$ 10,547
4	23,703	23,708	23,704	23,718
6	23,722	23,717	23,713	23,730
8	23,706	23,709	23,700	23,719

,

Table 3. Weight fractions (%) of PEG in production runs of all the simulated systems.

${\mathit{N}}_{\mathbf{PEG}}$	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36	72
$n=119,M_{\mathrm{r}}=5260$
8	8.98	8.96	8.94	8.89
12	12.89	12.85	12.83	12.77
16	16.47	16.44	16.41	16.33
$n=239,M_{\mathrm{r}}=$ 10,547
4	8.99	8.97	8.95	8.91
6	12.90	12.87	12.85	12.78
8	16.50	16.46	16.43	16.35

and

Table 4. Molalities of sodium chloride (mol/kg) in production runs of all the simulated systems.

${\mathit{N}}_{\mathbf{PEG}}$	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36	72
$n=119,M_{\mathrm{r}}=5260$
8	0.00000	0.0422	0.0844	0.1687
12	0.00000	0.0422	0.0843	0.1686
16	0.00000	0.0422	0.0844	0.1686
$n=239,M_{\mathrm{r}}=$ 10,547
4	0.00000	0.0421	0.0843	0.1685
6	0.00000	0.0421	0.0843	0.1684
8	0.00000	0.0421	0.0843	0.1685

. The number of water molecules, Table 2, was around 23,700, slightly varying between samples, mainly due to vapor deletion after the first equilibration phase. Concentrations of PEG, expressed as weight fractions, are given in Table 3. Finally, Table 4 lists molalities of sodium chloride, expressed as the number of moles of NaCl per 1 kg of water. Polymer concentrations were chosen in from the range where poly(ethylene oxide) of higher molecular weights (hundreds of kDa) spins [53,54], but correspondingly long chains could not be conveniently modeled by atomistic simulation. Thus, much lower molecular weight PEG was used, such as in our previous studies [33,34,35].

3.1.2. Visual Observation

Molecular dynamics simulations performed in the present study showed the key features of the electrospinning onset, namely the Taylor cone formation and jetting. The time evolution of the simulated systems under the field of 1.5 V/nm is illustrated by the sequence of sample frames along MD trajectories of two selected systems with eight 10.55-kDa chains: one without ions, see Figure 7, and the other with 72 NaCl ion pairs, see Figure 8. Both sequences show gradual deformation of the initially spherical droplet surface, progressing to the formation of the Taylor cone at ca. 30 ps followed by ejection of mass from the cone’s tip. Before a single, well-formed jet appears, several ‘tentacles’ of solvent molecules randomly protrude from the disrupted surface.

There are some notable differences between the systems with and without ions. At first glance, samples without ions appear to be more compact, while ions seem to increase surface roughness and to enhance the fragmentation of the liquid as the process advances. Sodium cations often enter the tentacles of the solvent described in the previous paragraph and cause them to break, escaping the droplet together with a number of water molecules around them. Figure 9 shows the comparison of a system without ions and a system with ions at the same time instant—escaping sodium ions accompanied by a number of solvent molecules are clearly visible. Another interesting observation can be made concerning the interaction of ions with PEG chains. It seems ions assist the unfolding of the chains and their transport to the jet. This might be caused by the coordination of PEG ether oxygens to Na${}^{+}$ ions which are dragged by the electric field, as suggested by Wang et al. [36].

3.1.3. Transfer of Mass and Jetting Onset Times

The total net mass transport in the direction of the field through the plane at $z=20$ nm by the given time instant in the production runs with the applied field strength of 1.5 V/nm is plotted in Figure 10 and Figure 11 for 5.26-kDa and 10.55-kDa chains, respectively. Figure 12 shows the mass transport of the system with eight 5.26-kDa PEG chains and different number of ions for different field strengths. Figure 13 and Figure 14 show the mass transport of PEG chains alone for systems with 5.26-kDa and 10.55-kDa PEG, respectively.

The onset of jetting was characterized by the production time at which mass transfer exceeded a certain threshold.

Table 5. Production times (ps) at which the total mass transfer through the plane at $z=20$ nm exceeded 2000 Da under 1.5 V/nm.

${\mathit{N}}_{\mathbf{PEG}}$	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36	72
$n=119,M_{\mathrm{r}}=5260$
8	42	46	48	42
12	50	42	50	48
16	50	50	50	38
$n=239,M_{\mathrm{r}}=$ 10,547
4	44	46	44	44
6	48	58	48	44
8	46	46	46	40

and Figure 15 show the times when total mass transported through $z=20$ nm exceeded 2000 Da, whereas

Table 6. Production times (ps) at which the mass transfer of PEG through the plane at $z=20$ nm exceeded 500 Da under 1.5 V/nm.

${\mathit{N}}_{\mathbf{PEG}}$	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36	72
$n=119,M_{\mathrm{r}}=5260$
8	>100	64	88	56
12	80	62	58	56
16	62	56	56	66
$n=239,M_{\mathrm{r}}=$ 10,547
4	>100	68	56	48
6	66	70	56	48
8	68	80	60	56

and Figure 16 show the times when mass of PEG alone transported through $z=20$ nm exceeded 500 Da The dependence of the times when total mass transfer exceeded 2000 Da on the intensity of the field is shown in

Table 7. Production times (ps) at which the total mass transfer through the plane at $z=20$ nm exceeded 2000 Da for different field intensities, E. Systems with eight 5.26-kDa chains were used.

E (V/nm)	${\mathit{N}}_{\mathbf{NaCl}}$
	0	18	36
1.0	82	84	84
1.5	42	46	48
2.0	30	30	30
2.5	24	24	24

and Figure 17 for the system with eight 5.26-kDa chains and different numbers of NaCl ion pairs.

An interesting observation is that the presence of ions had no significant impact on the total mass transfer and the corresponding onset times, while in most cases there was a discernible pattern of ions facilitating the transport of PEG chains alone, in accord with the visual observations presented above and the conclusion of Wang et al. [36]. This finding is further corroborated by Figure 18 and Figure 19 showing the velocity of PEG’s center of mass in the direction of the field, determined as a numerical derivative of the z-coordinate of the center of mass, see Equation (1). The accelerating effect of ions on PEG chains is clearly discernible, especially in the case of 10.55-kDa chains. As regards the effect of the electric field strength (Figure 17), the increase of mass transport and decrease of the corresponding onset times is not surprising and in agreement with findings at the macroscale [58].

3.1.4. Evolution of the Dipole Moment

Dielectric and ionic components of the total dipole moment under the field of 1.5 V/nm are plotted as functions of time in Figure 20 for 5.26-kDa PEG and Figure 21 for 10.55-kDa PEG. Curves for different field intensities are plotted in Figure 22 for selected systems. One can see that in all cases, dielectric polarization increased very steeply within the first 2 ps and then further moderately increased in contrast with ionic polarization, which is negligible in the beginning but takes off steeply in later stages of the process. The evolution of the dielectric dipole moment can be explained by fast dielectric saturation of molecules in the spherical droplet followed by further increase as the droplet deforms by polarization forces. The results on the dipole moment reinforce the findings that ions have little effect on mass transfer. It seems that at the molecular scale, ions play almost no role in the formation of the Taylor cone and the onset of jetting. In ion-containing systems, the same fundamental mechanisms apply as in systems with no ions at all, pointing to the clearly dielectric nature of the molecular-scale electrospinning onset.

3.1.5. Evaporation of Solvent

The rate of evaporation of solvent determined as an increase of the number of one-molecule water clusters through a chosen time period was evaluated for the periods 0–36 ps, 36–72 ps, and 0–72 ps of production time. Figure 23 shows the results for all 24 systems under 1.5 V/nm. The time 36 ps roughly corresponds to the instant when the Taylor cone has already fully developed and jetting starts under 1.5 V/nm, see Figure 7 and Figure 8. One can see that in the jetting phase (36–72 ps) and overall (0–72 ps) most simulations show the highest evaporation rate for the highest ionic concentration, but the trend is not always monotonous. Similarly to Wang et al. [36], the present authors suspect that the transfer of the kinetic energy from field-accelerated ions to the solvent is the main mechanism how the electrolyte contributes to the evaporation. In order to verify this assumption, the energy dissipated from ions was calculated for the given time periods and correlated with the evaporation rates, see bellow.

3.1.6. Energy Dissipation

The ionic energy dissipated from the beginning of production, $E_{\mathrm{dis}}(0,t)$, see Equation (2), as a function of time t, is shown in Figure 24 for selected systems, for which velocities along the trajectory were stored to allow for the evaluation of the kinetic energy.

In order to reveal a possible link between evaporation and energy dissipation, the evaporation rates in 0–36 ps, 36–72 ps, and 0–72 ps of production time were plotted against the energies dissipated through the same intervals, see Figure 25 and Figure 26. Linear regression was made and Pearson’s correlation coefficients, r, between the dissipated energy and the evaporation rate in the respective time intervals were also determined, separately for 5.26-kDa and 10.55-kDa PEG. For 5.26-kDa PEG, correlation coefficients were −0.3399, 0.9210, and 0.9125 for 0–36 ps, 36–72 ps, and 0–72 ps, respectively. For 10.55-kDa PEG, correlation coefficients were 0.3058, 0.7294, and 0.7307 for 0–36 ps, 36–72 ps, and 0–72 ps, respectively. The values of r in jetting phase and overall point to a clear correlation between the dissipated energy and the evaporation rate. Moreover, in all series of samples with the same number and length of PEG chains (same symbols in Figure 25 and Figure 26), there is always a non-decreasing trend in jetting phase and overall, unlike in the case of the dependency of the evaporation rate on the number of ion pairs, Figure 23.

3.2. Experimental

3.2.1. Preparation of Solutions

A total of 12 solutions were prepared using ultrapure water, PEG with number-average molecular weight 6000 Da, and NaCl. The composition of all solutions is summarized in

Table 8. Weight fractions (%) of PEG in the solutions used in experiments.

PEG Series	NaCl Series
	0	1	2	3
a	9.01	8.99	8.97	8.91
b	12.89	12.86	12.84	12.78
c	16.46	16.42	16.39	16.33

and

Table 9. Molalities of sodium chloride (moles of NaCl per kg of water) in solutions used in experiments.

PEG Series	NaCl Series
	0	1	2	3
a	0.0000	0.0429	0.0838	0.1694
b	0.0000	0.0427	0.0838	0.1693
c	0.0000	0.0428	0.0839	0.1695

and roughly corresponds to the composition of solutions used in molecular simulations, see the previous subsection.

Table 10. Conductivities of all samples ($\mathsf{\mu }$S/cm) at 21 ${}^{\circ }$C. Mean values calculated from 3 measurements are given along with standard deviations of the mean.

PEG Series	NaCl Series
	0	1	2	3
a	35 ± 1	2019 ± 54	3703 ± 25	6573 ± 31
b	37 ± 1	1845 ± 5	3343 ± 6	5923 ± 67
c	38 ± 1	1614 ± 4	2930 ± 10	5293 ± 21

lists the measured conductivities of the solutions.

3.2.2. Oscilloscopic Measurements

Table 11. The experimental instability onset times (ms) determined from oscilloscopic voltage–current delays. Mean values calculated from 10 measurements are given along with standard deviations of the mean. The air temperature was 24.0 ${}^{\circ }$C, and the relative humidity was 60%.

PEG Series	NaCl Series
	0	1	2	3
a	42.20 ± 2.16	41.92 ± 1.82	39.04 ± 2.87	35.60 ± 2.64
b	46.30 ± 2.24	39.60 ± 1.95	41.04 ± 3.26	39.96 ± 3.03
c	41.08 ± 2.30	40.16 ± 1.52	40.03 ± 1.79	42.20 ± 2.00

and Figure 27 summarize the results of voltage–current delay measurements used to characterize the onset of the instability of the meniscus. The measurements were done using setup I (see Table 1) at air temperature 24.0 ${}^{\circ }$C and relative humidity 60%. PEG series a and b show a weak tendency of the onset time to decrease with increasing ionic concentration, whereas the solutions with highest polymer concentration, PEG series c, do not show any trend.

Experimental onset times determined in this work show much weaker dependence on ionic concentration than expected for a process driven by the free charge density. For comparison, Yalcinkaya et al. [59,60] observed striking dependence of jet launching times on ionic concentration in needleless electrospinning of polyurethane and poly(ethylene oxide) solutions. Nevertheless, it should be pointed out that the instant of the first discharge measured in this work is not the same as the instant where jetting sets in. In order to elucidate what phenomena take place at the electrode, a high-speed camera coupled with oscilloscope was used.

3.2.3. High-Speed Photography

The progress of droplet deformation in setup I was recorded by a high-speed camera to complement the oscilloscopic data with visual image. Figure 28a shows the electrode current as a function of time for sample b0. The detail of the current onset phase is inserted. Figure 28b shows four photographs captured at the times marked in the detailed plot by vertical lines. Photographs show a well formed cone with an initially blunt apex, which sharpens when the first current peak is observed. No jet was visible during the onset phase, which could be caused by either the jet being too tiny to be captured in this setup or, more probably, by the absence of the jet. In the latter case, the authors believe that only the corona discharge [2] from a sharp tip of the Taylor cone was responsible for the current signal, but more data are needed to confirm this scenario. A visible jet was captured in later stages of the process with sample a0, see Figure 29. One can clearly see the ejection of a thin liquid jet from the tip of the Taylor cone and its eventual decay into visible droplets. From the time axis of Figure 29, one can get an idea how the launching of jets is delayed with respect to the first discharge, which occurred at ca. 35 ms for this particular measurement. The lateral position of the cone and the jet is caused by inhomogeneous electric field, which strengthens near the edge of the electrode.

In setup II, with the small-diameter electrode positioned horizontally, the first discharge was always accompanied by the detachment of a large portion of liquid in the form of a ‘spindle’ and large droplets. See Figure 30 for an example of the spindle observed with solution a0. Jetting was more prevalent than in setup I. A thin jet is also observable in Figure 30 at the tip of the spindle. Figure 31 and Figure 32 show cone-jet formations in later stages of the process. In Figure 31 the jet eventually transforms into visible droplets, whereas in Figure 32 no droplets are apparent. Noticeable splaying of the jet in Figure 32 might be caused by whipping instability prior to disintegration to aerosol.

3.3. Remarks on Scales

There seems to be a fundamental difference between macroscopic theories of electrospinning and the atomistic picture presented in this work. It is well established in macroscopic continuum models that electrospinning, especially in the Taylor cone formation stage, is governed by mobile (ionic) charges, and that even the fluids with the slightest concentrations of ions, unavoidable in any liquid, can be considered conducting in this regard [5]. A hypothetical pure dielectric fluid is also predicted to form a conical meniscus but not able to jet, though [61]. The above macroscopic findings seem to be in contradiction to what is observed in atomistic modeling, where ions play minor role and do not significantly affect the nature of the field-induced transformation of a dielectric liquid. An interesting analogy is offered by experiments with ferrofluids, which do jet in magnetic fields [62]. Ferrofluids can be viewed as purely dipolar in principle due to nonexistence of magnetic monopoles.

The above discussed disparity deserves some explanation. First we have to realize that in simulations we work at time scales many orders of magnitude smaller than in experiments. Simulation times in tens of picoseconds put the considerations on charge relaxation to completely different perspective. The length scales are also completely different from regular electrospinning. Ions do not have enough time to create a charged interface and much possibly enough space to rise proper capillary waves against surface tension. That might shift the jetting process to the one governed by polarization forces. Hartman et al. [63] calculated with their model that although the surface charge of jetting liquid predominantly comprised the free charge in most parts of the surface of the cone and the jet, the ratio of polarization charge sharply increased at the point of the cone–jet junction, suggesting at least some role of polarization forces in this region. This lends support to the interpretation of the content of the simulation cell as a very apex of the Taylor cone where the jet starts, in accord with the simulation methodology proposed by Jirsák et al. [31] and followed herein. Contracted time and length scales are also coupled with much higher electric fields needed to observe the process. The intensities of electric fields used in simulations of electrospinning and other related phenomena are usually much higher than the average values attainable in experiments [31]. We note in passing that there are reports on related phenomena which approach the typical size of the simulation, e.g., the water bridge between an AFM tip and a sample [64].

3.4. Limitations of the Study

The present study has several limitations that should be understood for correct interpretation of the results. Some results of MD simulations lack statistical power due to limited number of simulated trajectories. The authors preferred using computing resources to cover a broader concentration range rather than to obtain several independent runs for a limited number of systems. Only nonpolarizable models were used for all components, which could affect dynamic properties. It is known that nonpolarizable models of ions and water do not correctly reproduce water diffusivity in electrolyte solutions [65]. Similarly to all molecular simulation studies in electric fields, it is not clear how to exactly relate the field applied in simulation to the macroscopic one.

As regards laboratory experiments, it should be pointed out that they play only a complementary role in the present study and were mainly motivated by authors’ curiosity to learn how the solutions used in simulations would present themselves at macroscale. The solutions were expected to spray rather then spin, which showed to be the case. The most important information the authors wanted to extract from the experiments was the dependence of the characteristic time of the Taylor cone build-up on composition. In future studies, a search should be made for the set of conditions where jetting in simulations is more closely related to real fiber-producing experiments. The technique of measurements could also be improved. In the present study, the volume of droplets and their exact position with respect to the edge of the electrode were not precisely controlled for, which introduced additional scatter to the measured onset times, as the field intensity varies very steeply near the edges.

4. Conclusions

Molecular dynamics simulations of a nanosized droplet of polymer solution subjected to strong electric fields were performed. Formation of the Taylor cone and subsequent jetting accompanied by rapid evaporation of solvent was observed. The results indicate that at molecular scale, the formation of the cone and the onset of jetting is governed by dielectric forces, in apparent contrast with macroscopic continuum models. Ions do not significantly affected the nature of the process but facilitated the transport of polymer chains into the jet and enhanced the evaporation of solvent.

Experiments were performed with similar composition of solutions as in the simulations. A weaker-than-expected effect of ions on the time of formation of the Taylor cone was found—onset times slightly decreased with increasing ionic concentration, with the exception of the highest polymer concentration where no trend was discernible. High-speed photographs show that the solution did not actually spin, which was expected. The observations are consistent with electrospray in different modes (corona discharge, cone–jet, spindle). There is not enough data to assess the effect of ions on the morphology of the process.

Apparent discord between governing mechanisms, i.e., conductive in macroscale models vs. dielectric in molecular models might be upsetting, but creates an opportunity to grasp better understanding of the underlying physics in its full complexity. Simple considerations on the time and length scales involved in experiments and in simulations might bring a relief that there is nothing fundamentally wrong with either, but only after attempting to bridge the vast space between nano- and macroscale we may truly learn something new. That is the path the present authors strive to follow in future research.