Electric fields can control the transport of water in carbon nanotubes

The properties of water confined inside nanotubes are of considerable scientific and technological interest. We use molecular dynamics to investigate the structure and average orientation of water flowing within a carbon nanotube. We find that water exhibits biaxial paranematic liquid crystal ordering both within the nanotube and close to its ends. This preferred molecular ordering is enhanced when an axial electric field is applied, affecting the water flow rate through the nanotube. A spatially patterned electric field can minimize nanotube entrance effects and significantly increase the flow rate.


Introduction
Recent advances in materials science have meant that the study of water within carbon nanotubes (CNTs) is now an important research topic [1]. In particular, computational and experimental studies of water flow inside aligned CNT membranes have shown the potential for creating filtration membranes with very promising capabilities [2].
Very narrow CNTs, of approximately 0.8 nm diameter, have attracted particularly strong interest as they resemble biological water channels such as aquaporins, proton pumps and protein cavities in lysozyme or Alzheimer's Aβ amyloids. The filling process, and the structure and orientation of water inside this type of  biaxiality, it is possible for both uniaxial and biaxial molecules to exhibit biaxial ordering [17]. A water molecule is molecularly biaxial, i.e. the molecular orientation must be defined by an orthogonal set of molecular directions. It is perhaps the simplest biaxial molecule that naturally occurs, so a biaxial ordering might be expected to be possible. In previous studies, the structure and average orientation of water close to solid surfaces or inside CNTs have usually been presented in terms of the density, radial distributions, hydrogen bonding, dipole moments, entropy, enthalpy and free energy [18][19][20][21][22]. These measures provide an adequate description of the spatial structuring of the water, and the presence and position of layers. The angle between the direction l m (the dipole moment vector) and an axis (such as the flow direction or an electric field direction) defines the dipole moment, the distribution of which can provide an average orientation and a measure of directional ordering. In addition, ordering of confined water in CNTs has been described in the past in terms of entropy, where decreased rotational entropy highlighted the increased ordering of water [19].
However, for a more complete description of the ordering of biaxial molecules, three orthogonal molecular directions n m , l m and m m must be defined, such as those for the water molecule shown in figure 1. The direction n m corresponds, by convention in the liquid crystal field, to the molecular long axis; l m is along the molecular short axis (which coincides with the molecular dipole direction); and m m is orthogonal to n m and l m (i.e. m m = n m × l m ). With these directions, the symmetric, traceless, second-rank tensors can be constructed: where I is the second-rank identity tensor. The average orientation of an ensemble of these biaxial molecules can then be represented by the two order tensors Q and B, which are averages of their microscopic counterparts q and b. We consider both the radial (measured from the nanotube axis) and axial (along the nanotube axis) dependence of the order tensors, i.e. Q(r) and Q(z) (we find no statistically significant dependence of the order tensors on the polar angle θ ). Ensemble averages are constructed over bins in the streamwise (axial) direction or the radial direction from (2.4) where N is the number of water molecules inside the axial or radial bin andτ is the number of MD time steps over which the system is sampled. When the water density ρ is small, the statistics of Q and B will be unreliable, so only measurements for ρ > 100 kg m −3 are considered here. We will show that these order tensors provide information beyond a standard dipole moment calculation, the information for which is in fact contained within one of the eigenvalue/eigenvector pairs of the tensor B.
The measured values of Q and B are interpreted through an investigation of their eigenvalues and eigenvectors. The tensors are symmetric and traceless so that two of the total three eigenvalues are independent, and the set of eigenvectors forms an orthogonal triad, for each tensor. It is the eigenvectors that describe the principal axes of the averaged molecular axes, and the eigenvalues measure the spread of molecular deviation from the eigenvector: Q provides information about the long molecular axis, and B describes the effects of molecular biaxiality on the ensemble average properties. Both tensors are zero in the isotropic phase (i.e. bulk water), but any loss of isotropy will lead to non-zero eigenvalues and a sensitivity to the direction of the average orientation of the molecular directions.
If the molecular direction n m forms a uniaxial arrangement then Q will have two equal eigenvalues and a third distinct one; similarly for B. In uniaxial cases, the eigenvector corresponding to the distinct eigenvalue indicates the orientation of the director, revealing a radial, angular or axial uniaxial nematic phase. If the distinct eigenvalue is negative then the molecular axis is, on average, oriented randomly in a plane with a director oriented orthogonal to this plane.
On the other hand, biaxiality is present if all three eigenvalues of either Q or B are different, with the former case indicating that the molecular direction n m has a biaxial configuration and the latter case revealing a biaxial nematic phase due to the molecular shape. As the bulk phase of water is isotropic, if there is any form of anisotropic order, uniaxial or biaxial, it has been induced by the confining CNT or an applied electric field, or both. Such a situation is termed 'paranematic ordering'.
Convenient, and commonly used, measures of the presence of anisotropy and biaxiality are through the parameters [23,24]:

Simulation method
We carry out MD simulations of water flow through CNT membranes with an applied pressure drop and three different applied electric fields, as illustrated in figure 2. In all our cases, we simulate a 12.5 nm long rigid CNT (which is a semiconductor) that connects two water reservoirs at the ends. Two rigid graphene sheets of dimensions 10.5 × 10.2 nm, and each with a pore equal to the CNT diameter, are placed at each end of the tube to act as separating membrane walls.
The chosen CNT in this study is of diameter D = 2.034 nm, chirality (15,14) and 0.142 nm carboncarbon bond length. The CNT is long enough to produce a fully developed water flow inside it within a reasonable simulation time. Our MD simulations are performed using the mdFOAM solver [25,26], which has been developed within the open source OpenFOAM software (CFD Direct Ltd, http://www.openfoam. org). This has been used previously to investigate various other micro/nano flow problems, including flows of water through CNTs [8,[27][28][29]. All molecules in a MD simulation evolve in rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: and where i = (1, 2, . . . , N) is the index of a molecule in a system of N molecules, and m i , r i , v i , f i are the molecule's mass, position, velocity and total experienced force, respectively. The total force on a molecule is computed from a sum of pair intermolecular forces with neighbouring molecules by The pairwise potential U(r ij ) contains a Lennard-Jones (LJ) potential and Coulomb interactions for charged sites, as indicated by the first and second terms on the right-hand side of the following equation: where r ij = |r i − r j | is the separation of two arbitrary atomic sites (i, j) of species (a, b), σ ab and ab are the LJ characteristic diameter and well depth, respectively, q a is the electric charge of site a, and ε 0 is the permittivity of vacuum. Equations (3.1) and (3.2) are integrated numerically using the Verlet algorithm with a discretized time step t = 2 fs, while the pair force calculation step in equation (3.3) is computationally optimized by using double-loop savings, a cut-off distance r cut , and the cell-list algorithm. Both LJ and Coulomb interactions are truncated at a cut-off radius of r cut = 1.0 nm [30], while the Coulomb interactions are also shifted [31,32].
We use the TIP4P/2005 water model [33], which consists of four atomic sites: one oxygen O, two hydrogen H, each with charge q H = 0.5564e, and one massless M located 0.1546 Å in the direction of l m (figure 1) with charge q H = −1.1128e. The oxygen sites have zero charge and as a consequence they only interact through the LJ potential, with the parameters σ OO = 3.1589 Å and OO = 1.2868 × 10 −21 J [33].
The carbon-water potential consists only of a carbon-oxygen (C-O) LJ interaction, given by the parameters σ CO = 3.19 Å and CO = 7.09302 × 10 −22 J that reproduce the macroscopic contact angle of a water droplet on a graphitic surface [8]. We use Hamilton's quaternions in order to keep the fixed geometry of the water molecules, with an O-H bond distance of 9.572 Å and an H-O-H angle of 104.52 • . The solid atoms in the CNT and the graphitic sheets are kept fixed, which is a common simplification in MD simulations [6,34].
All our MD simulations are in an NVT ensemble. Periodic boundary conditions are applied in all directions of the domain, and a Berendsen thermostat is applied only in the two reservoirs to maintain an overall water temperature of 298 K. Every case is initially equilibrated for at least 4 ns of problem time. During the first 800 ps of the equilibration run, density is also controlled inside the reservoirs (indicated by blue boxed regions in figure 2) using the FADE insertion/deletion protocol [35]. The target densities can be calculated according to the compressibility equation: where ρ 0 = 997.5 kg m −3 is the density of water under atmospheric pressure at 298 K, E bm = 2.15 × 10 9 N m −2 is the bulk modulus fluid elasticity of water, P 1 = 100 MPa and P 0 = 0.1 MPa are the high and low (atmospheric) pressures in the corresponding reservoirs, respectively. The calculated density from equation (3.5) for the high-pressure reservoir is 1046 kg m −3 . After equilibration, each simulation contains approximately 33 000 water molecules. Production runs are then performed for a minimum duration of 4 ns, during which all measurements are sampled and averaged in order to achieve better statistics. This required at least 110 computational hours on a graphical processor unit (GPU). In all the simulations carried out in this paper, a hydrostatic pressure drop is imposed across the two reservoirs. This is enabled by applying a force to all water molecules that reside at any point in time during the simulation in a slab (of thickness z = 3 nm) centred at the outer part of the reservoir (z = 0), with the forcing Gaussianly distributed [36], i.e.
where σ s = z/8 is the chosen standard deviation of the Gaussian, ρ n is the fluid number density, F z (z) = 0 outside the forcing region and P = 100 MPa is the chosen pressure drop. This pressure difference, although high by experimental standards, is common practice in MD in order to increase the signal of the flow rate above the thermal noise of the molecules, and as a consequence reduce the computational cost of the simulations [30,37,38]. The applied pressure difference creates a net mass flow rate ṁ (kg s −1 ) along the nanotube, which is measured at a plane in the middle of the tube. The mass flow rate is calculated as a sum over the net molecules that cross the flux measurement plane in the z-direction: where the first summation is over the time period t → t +τ t, and the second summation is over water molecules that cross the plane in one time step. The majority of measurables presented in the Results section ( §4) are taken inside the CNT, between z = 5.7 nm and z = 15.8 nm.
In the three cases indicated in figure 2, in addition to the imposed hydrostatic pressure drop, we also apply an external homogeneous static electric field E along the z-axis of the system and l Q [1] l Q [2] l Q [3] r l B [1] l B [2] l B [3] [3] l B [1] l B [3] (c) Figure 3. (a) Water density profile and eigenvalues λ Q of the order tensor Q, and (b) water density profile and eigenvalues λ B of the order tensor B, as a function of radial distance from the centre of the CNT axis, averaged over the length of the CNT for the zero-field situation. (c) Maximum and minimum eigenvalues λ Q,B of the order tensors when the electric field is applied over the whole system (Case A). The standard error of the mean (s.e.m.) with 95% confidence interval (CI) is smaller than the symbol. (Online version in colour.) localized to specific regions of the domain. We note here that the imposed electric fields are somewhat artificial, as the field lines in reality can be curved by the differences in permittivities in the liquid (if not pure water) and between the liquid and the membrane. Taking this into account, our results could be slightly affected by more realistic fields.
The electric field is implemented as an additional force F i = q i · E to atoms within the water molecule, where q i is the partial charge on atom i [39,40]. As water molecules are neutral but polarized, the effect of the electric field on each molecule is a torque. It has previously been shown that this is a valid MD representation of the influence of an electromotive force exerted by a voltage difference between two electrodes [41,42]. This method has also been used to impose a potential difference across a membrane in simulations of electroporation of lipid bilayers [43].

Results
We consider the anisotropy and biaxiality of water transported along a CNT (i.e. in the z-direction) by a pressure difference, and the subsequent changes in ordering when an axial electric field is applied. The electric field strengths we examine range from 0.01 to 1.0 V nm −1 , which are typical of those produced by charged electrodes, ion channels, ionic biomolecules or assemblies, and electroporation [44,45]. In biological membranes during the process of electroporation, for example, the typical electric field range is 0.3-0.5 V nm −1 [46,47].
When an electric field in the range 0.01 to 1.0 V nm −1 is applied to water in a nanotube, the density and hydrogen bonding profiles are all almost identical to a reference zero-field case. A number of previous MD studies have used much stronger electric fields (greater than 1.5 V nm −1 ) to produce in the water a structural ordering similar to that of ice [20,22,48]. Other  x Q [2] x Q [3] x B [1] x B [2] x B [3]   recent studies have shown applied electric fields inducing or promoting the liquid-to-solid phase transition [40,49]. In all the cases we present here, however, the water remains in a liquid phase.

(a) Water confined in nanotubes
We first perform pressure-driven simulations without an applied electric field: this is our zerofield reference case. In the reservoirs, and within the bulk region inside the nanotube, our results show that water is isotropic (α Q ≈ 0, α B ≈ 0). However, we find that water close to the CNT walls is in a biaxial paranematic state. Figure 3a shows this weak uniaxial state for the molecular direction n m , because Q has two equal eigenvalues and one distinct, at a radial distance of approximately 0.65 nm from the CNT axis, which is where the density is at a maximum. Figure 3b shows that, at this same location, the B tensor exhibits a degree of biaxiality, as all the eigenvalues are different.
If we now apply an electric field over the entire system, anisotropy rises and water molecules start aligning their dipole moment vector along the same direction as the applied electric field.  x B [1] x B [2] x B [3] (e) ( f ) Figure 5.
(a-f ) The z, r, θ components of the eigenvectors ξ for all the eigenvalues of the order tensors Q and B for Case A, with an applied electric field across the entire system of (0, 0, 0.959) V nm −1 . The radial density profile is also plotted, in blue, in order to highlight water layering. Eigenvector components are plotted only when the corresponding eigenvalue is greater than 0.02, in order to highlight areas with significant anisotropy and to reduce noise. Note that eigenvectors are defining symmetry axes, rather than directions along axes. The eigenvector component values of +1 and −1 are equivalent and any abrupt change between these extreme values is only an artefact of the use of vector components. (Online version in colour.) For the largest applied field (E = (0, 0, 0.959) V nm −1 ), λ Q [1] still has the highest absolute value among the Q eigenvalues, while λ B [3] has the highest absolute value of the B eigenvalues, as figure 3c shows. An investigation of the eigenvectors of the tensors for the zero-field case, shown in figure 4a,b, indicates that the preferred orientation of l m is radial to the nanotube axis. In detail, for both tensors, the first eigenvalue (λ Q [1] and λ B [1] ) has the largest absolute value, meaning that eigenvectors ξ Q [1] and ξ B [1] will reveal the preferred average direction of the water molecules inside the CNT. In both eigenvectors, the radial component (black diamonds on figure 4) is the dominant one within the high-density and high-anisotropy region near the wall. At first, this may seem contradictory, as Q and B are orthogonal by construction. But as both of the largest absolute eigenvalues (λ Q [1] and λ B [1] ) are negative ( figure 4a,b), then the molecular axis is, on average, oriented randomly in a plane, with the eigenvector oriented orthogonal to the plane.  It can therefore be concluded that the average preferred orientation of the dipole moment vector l m is radial to the nanotube axis.
When an electric field is applied, the eigenvectors ξ Q [1] and ξ B [3] , shown in figure 5, clearly show that the average preferred orientation of l m has changed from radial to the nanotube axis to axial along the CNT, in the direction of both the flow and the applied electric field.
In order to further analyse the structure of water molecules inside the CNT without an external electric field, radial profiles of density and the average number of hydrogen bonds (n HB ) are presented in figure 6. Towards the centre of the CNT both the density and hydrogen bonding have values close to the bulk values. However, in the high-density water layer formed close  to the CNT walls, the number of hydrogen bonds is significantly reduced, as expected, due to the presence of the solid wall. This highly dense and oriented water layer promotes water flow inside the nanotubes, essentially acting as a lubrication film between the bulk water phase and the solid walls.

(b) The effect of patterned electric field configurations
When an electric field is applied to the system, its orienting effect on the molecular dipole alters the fluid anisotropy and biaxiality. However, as reported in previous publications [20,39,48], the structural nature and liquid phase of water under these conditions is not altered, so it is impossible to reveal the average orientational change of the water molecules from density and hydrogen bonding measurements alone.
We consider three regions over which we apply the electric field. In the first situation, Case A, all molecules in the system experience the electric field and we see, in figure 7a-c, that anisotropy (α Q,B ) increases with increasing field strength. The biaxiality of Q also increases (figure 7b), but the biaxiality of B initially increases but then decreases as the field strength increases (figure 7d). This decrease in biaxiality is due to the aligning effects of the applied field (a forcing that has uniaxial symmetry), which tends to promote uniaxial ordering of the molecular dipole along the z-axis, in competition with the zero-field alignment of l m , which is in the radial direction. In figure 8, we see that in Case A the mass fluxṁ E , normalized by the zero-field fluxṁ 0 , does not appreciably change. Although anisotropy is increasing significantly as we increase the field strength, the mass flow rate remains almost constant for E < 0.7 V nm −1 . For higher electric fields up to E = 0.959 V nm −1 , the predicted flow rate increases only by approximately 10% compared with the zero-field case.
In Case B, where the electric field is applied over the CNT region only (i.e. from z = 4.7 nm to z = 16.7 nm), the anisotropy and biaxiality parameters are similar to Case A. However, in this case, figure 8 shows a decrease in the mass flow rate with increasing field strength-by approximately 30% at the highest field strength values. The changes in order parameter within the CNT are causing an increased dissipation in the system, particularly around the nanotube entrance.
In Case C, the electric field is only applied in the entrance region of the CNT (from z = 3.2 nm to z = 4.7 nm). The mean values of anisotropy remain relatively low, as the orienting field only affects a limited region. Biaxiality does seem to occur but the error bars are large because the anisotropy (which occurs in the denominator of β Q,B 2 ) is small. To highlight the differences between Cases   electric field E = (0, 0, 0.384) V nm −1 in both cases. The anisotropy rise at the entrance of the CNT in Case C, in contrast with the high anisotropy along the whole of the CNT in Case A, seems to significantly affect the flow rate. Figure 8 shows that this localized application of an electric field leads to a linear increase of the flow rate with field strength. At the highest field strength studied, we measured flow rates of more than two-and-a-half times the zero-field flow rate, even though the cross-membrane pressure difference remains the same.
We have already shown in figure 7 that the orienting effect on the water molecular dipole of an electric field applied anywhere in the system will alter the anisotropy and biaxiality. However, the density structuring and liquid phase of water under these conditions are not altered, and so it would be impossible to uncover this change in the average orientation of the water molecules by examining the density and hydrogen bonding measurements alone. In fact, for  the range of electric field strengths applied here, the shapes of the radial density and hydrogen bonding profiles are nearly identical to those in the zero-field case, as a few characteristic field magnitudes show in figure 10a,b. Although figure 10a shows that the peak and trough values of the density profile change as a stronger electric field is applied, the average density inside the CNT is altered by no more than 2%. Similarly, for the hydrogen bonding measurements in figure 10b, the maximum difference between the high-field results and the zero-field results is even smaller than the maximum density difference (and in fact both measures are close to the level of statistical error).
To understand the orienting effects of the spatially patterned fields in each of the cases, we consider the variations in pressure, velocity, anisotropy and biaxiality along the axial direction z. Figures 11 and 12 show plots of these profiles through the system (i.e. from the inlet reservoir, through the CNT, to the outlet reservoir) for the zero-field case, as well as Cases B and C, which were shown to reduce and increase flow rates, respectively.
These results indicate that two important effects are determining the behaviour. Previous work has shown that increased isotropic order can increase the viscosity, while increased biaxial order may either increase or decrease the viscosity depending on the flow profile present [50]. In a previous publication [8], we have measured the viscosity of water inside nanotubes of different materials and we found that close to the walls it is lower than in the bulk. These lower values of viscosity can now be attributed to the biaxial order of the water molecules near the walls. However, reorganization of the molecular ordering around the entrance of the CNT-as molecules obtain biaxial ordering before entering the CNT-is also important: it influences the elastic contribution to pressure, as seen in continuum approaches such as the Ericksen-Leslie theory [13]. While viscosity and elasticity are properties of a continuum, rather than the molecular system considered here, their microscopic origins play an important role. Such ordering effects are known to be affected by the application of an electric field [51], as we also see in this work. Figure 11 shows that the CNT naturally transports water with minimal dissipation, although entrance losses are quite significant. An electric field over the nanotube alone (Case B) not only increases the anisotropy and biaxiality of the water (figure 12) but also promotes uniaxial ordering of the molecular dipole along the z-axis, which changes the viscosity and dissipation, and leads to a reduced flow rate.  Figure 12 shows that, in Case C, the ordering within the nanotube is unchanged from the zero-field case, but anisotropy and biaxiality are increased close to the entrance of the CNT. This reduces the high deformation experienced by the molecular fluid as it enters the CNT, and decreases the pressure difference between the entrance region and the nanotube. This significant reduction in entrance pressure loss then enables the flow to be driven in the nanotube with a higher pressure ( figure 11). At the exit of the nanotube, expansion to the reservoir leads to a sudden decrease in pressure and velocity. Although this expansion is more severe in Case C, the effect on the flow velocity is minimal compared with the gains in velocity due to the reduced entrance losses.
These results suggest that an applied electric field produces competing effects on the mass flow rate: reducing flow when the field affects molecules within the CNT, but increasing flow rates when the field aligns molecules close to the inlet of the CNT. It therefore appears that an electric field can be used to 'pre-order' water molecules to form a paranematic biaxial phase in a region close to the CNT entrance. This consequently reduces entrance losses and produces much higher flow rates for the same applied pressure difference.

(c) Electric field and flow enhancement limitations
Finally, we investigate the effect of even stronger fields in Case C simulations, and discuss the maximum flow rate enhancement achievable by using fields of very high magnitude. These fields are much stronger than those commonly observed in biological processes or in electromechanical devices. It was mentioned above that extreme fields (greater than 1.5 V nm −1 ) do change the properties of water from liquid to solid (crystalline), and can induce or enhance the liquid-to-solid transition in low-temperature environments. Figure 13 shows our predicted maximum flow rate increase that can be achieved using a static electric field applied close to the CNT entrance (Case C). For fields greater than 2 V nm −1 , the water state change (i.e. solidification) hinders any gains from the pre-ordering effect we have  The water phase change can be verified through the density and hydrogen bonding measurements presented in figure 14. There is a sudden density decrease in the region where the electric field is applied, while hydrogen bonding increases. For electric fields greater than 1 V nm −1 , the density decrease is more significant, to approximately 950 kg m −3 -the density of ice and not liquid water. In addition, the average number of hydrogen bonds per water molecule increases from 3.6 in the no-field region to 3.7. For fields stronger than 3 V nm −1 there is a further increase, indicating the complete solidification of water in that region. In combination with the flow rate results in figure 13, it can be concluded that the optimum flow conditions are achieved when the fluid density is low but the number of hydrogen bonds has not increased significantly.
Anisotropy and biaxiality also increase in the region where the electric field is applied, as figure 15 shows. Figure 15a,c shows that, for field magnitudes up to 3 V nm −1 , the maximum value of anisotropy is measured inside the entrance of the CNT, while for higher electric fields the anisotropy is higher outside the CNT within the electric field region. The shape of the anisotropy curve is altered and the affected region becomes wider when the field is stronger than 3 V nm −1 . Although biaxiality is highly affected (figure 15b,d), the increase is not monotonic with the magnitude of the applied field. The value of parameter β 2 initially decreases as the electric field increases from 0.959 V nm −1 to 3 V nm −1 . A further increase in the electric field strength then significantly increases β 2 .

Summary
We have presented a method that offers a more complete description of ensemble order at the molecular scale, for application to nano flow problems. Using the order tensors of this method, an anisotropic/biaxial ordering of water inside and close to the entrance of a CNT was revealed in MD simulations. We predicted that spatially patterned axial electric fields could be used to either enhance or reduce the pressure-driven flow of water along nanotubes.
Our results indicate an approximately linear growth of flow rate with electric field strength up to 2 V nm −1 . At this field magnitude an increase of over 300% in the flow rate, with respect to the zero-field case, was observed. There are, of course, physical limits to the magnitude of the electric field that can be applied, and a further increase induced a phase change in the fluid that hinders or even prohibits flow through the CNT. The value of this maximum field strength may well be affected by some factors we have not studied here, such as the geometry, chirality or composition of the nanotube. These aspects could be the subject of useful future research. In addition, other methods for inducing advantageous molecular alignment and ordering might be