Surface Ordering in Molecular Clusters by Interactions of Well-Defined Parity Multipoles.

The spontaneous alignment of dipoles in patterns tangential or transversal to the interface of molecular clusters is studied by molecular dynamics simulations throughout the entire aggregation process. Tangential ordering (TANGO) is found to rely on dispersion forces driven by dipolar fluctuations. Transversal ordering (TRANSO) results from the interplay of two conditions: the broken translational invariance at the cluster interface and the interactions of multipoles of opposite parity. In this case, the orientational order parameter at the interface follows a universal scaling law behavior. In the limiting case of strongly interacting particles, a sharp structural transition from a disordered into an ordered state is observed at a critical distance inside the cluster and at a critical magnitude of the dipole moment.

Introduction

Molecular clusters have attracted considerable interest in the past years[1] as important prototypes of the so-called small thermodynamic systems[2] which, under special conditions, can undergo phase transitions.[3] Relying on their ability to well reproduce ad incipio thermodynamic processes, they have been proposed as generic models for soft matter[4] and precursors for functional materials.[5]

Clusters of polar molecules can well reproduce the general process of aggregation and other more complex processes such as nucleation, solvation, and interface formation. The establishment of order at an interface is perhaps the most basic effect in determining the macroscopic properties of molecular clusters, such as polarization, surface entropy, and surface tension.

The present work aims at understanding the basic physics driving surface ordering and especially the tendency of polar molecules to align with respect to an interface. We perform large-scale molecular dynamics (MD) simulations of two-dimensional systems of biatomic molecules whose atoms are coupled by the Lennard-Jones (LJ) and the Coulomb potential in multipole expansion up to the second order. The molecules are neutral, so, at the lowest order, only the dipole and quadrupole moments enter the electrostatic force. Starting from an arbitrary initial configuration, molecules assemble in clusters. Parameters in the forces can be tailored so that the system evolves in a single cluster. The possibility of tuning the force field and the substantial freedom in choosing the initial condition, compatibly with the equilibrium temperature, open up an exciting research area where, in principle, the self-assembly of molecules can be engineered in such a way that the final nanocrystal has the desired physical properties. The choice of this very force field is motivated by the sound evidence that transversal ordering (TRANSO), that is, dipoles aligned transversally to the interface, results from the interactions of multipoles of different parity.[6−8] In this respect, our work differs from the vast majority of previous studies restricted to van der Waals and dipolar interactions only.[9] From a formal standpoint, it can rigorously be proven[10] that a volume distribution of dipoles (resp. quadrupoles) generates a surface density of charge (resp. dipoles) at the boundary, giving a sound basis to our assumptions.

The other possible pattern of orientational ordering exhibits molecules oriented tangentially to the interface and organized into end-to-end or end-to-tail shaped chains, which is the other major object of this study. We do not limit the analysis to the external interface of the system but attempt to explore the surface order well inside the system, at distances where the effect of the interface fades out but a bulk regime is not yet entered in a limiting sense. As a result, very specific length scales appear making molecular clusters very peculiar.

From an experimental standpoint, several techniques are able to resolve orientational order at the interface. Optical second harmonic (SH), sum-frequency generation (SH),[11] and X-ray reflectivity[12] are very established methods and have extensively been employed.

Model

We restrict ourselves to two-dimensional systems of N biatomic, rigid, neutral molecules. As such, they are modeled as a rigid body with an inertial moment arising from the two identical atomic masses in a dumbbell-like arrangement and with an intrinsic (point) electric dipole and quadrupole.

The kinematics of a molecule is thus described by the position of its center of mass (cm) and by a unit vector e giving the orientation. Even if we assume the two atoms to have identical physical properties in the model, it is useful for our discussion to distinguish them on the basis of the molecular orientation. We therefore label them as type “I” and type “II” if e points from “II” toward “I”.

Intermolecular interactions are modeled by the LJ and Coulomb potential. The latter is given in a multipolar expansion up to the second order. To fully understand how the different nature of the interactions may affect the results, the code developed is capable to selectively include or exclude each of the electrostatic forces in a way that will be clarified later.

The electrostatic energy of two molecules with the center of mass located in r and r is, at the second-order expansion, given bywhere the terms on the right-hand side are the dipole–dipole, dipole–quadrupole, quadrupole–dipole, and quadrupole–quadrupole potential energies, respectively, given bywhere r = r – r, r = |r|, θ (resp. θ) is the angle between r and the h-coordinate (resp. k-coordinate) direction, and γ is the angle between r and r. We emphasize that the following treatment is not restricted to any specific potential as long as it can be expanded in a series of functions of well-defined parity. The electrostatic force on molecule k due to h is given by[13,14]and the corresponding torque is given by[13,14]

The LJ potential energy of the two molecules in r and r is the superposition of the atom–atom contributions given bywhere r and r are the positions of atom i of molecule h and of atom j of molecule k, respectively. The two empirical parameters ϵ and σ that set the minimum energy and the range of the interaction, respectively, are given by[15,16]

We assume ϵI = ϵII and, if not otherwise specified, σI = σII.

The resultant force on molecule k is given bywhere FLJ(r) = −∇|VLJ. The resultant torque is given by

The general equation of motion of a rigid body readswhere F and τ are the resultant force and torque, respectively, is the inertia tensor relative to the center of mass, m is the mass of the molecule, and ω and α are its angular velocity and acceleration, respectively. If we are considering a planar motion, rotations can only occur around an axis orthogonal to the plane; consequently, the above equations simplify considerably. Let us assume, for instance, the motion to be confined in the xy plane, then only the z-component of ω and α is nonzero and eq for molecule k can be written aswhere is now the moment of inertia about an axis orthogonal to the plane and through the center of mass of the molecule. For linear molecules, the equations of motion can be further simplified because one of the components of the principal inertia moment vanishes and the other two are equal, specifically, in our configuration, I = I ≡ I and I = 0. Although the equations of the translational motion of the center of mass (eq ) can be straightforwardly integrated by conventional finite difference method,[17] the rotational ones require some additional manipulation to transform the second-order differential equations into a first-order system, given by[18]where e is the orientation vector of molecule k, g⊥ is the component of the torque (in the cm frame of reference) orthogonal to e, and λ is a Lagrange multiplier which constrains the dynamics to preserve the bond length. Equation can numerically be solved, for instance, by the Fincham algorithm.[19] The dynamics of the system is, finally, completely defined by the set of 6N equations, eqs and 15, with the initial conditions

The initial positions in eq are assigned by placing the cm of each molecule at the site of a two-dimensional lattice of size d ± δd, where δd is a random number different at each site and d is an input parameter which gives the initial average distance of the molecules, fixing, in this way, the strength of the initial force field. The initial orientations in eq are assigned by further rotating each molecule around an axis through the cm by a random angle. The center of mass initial velocity field is approximated by a Maxwellian distribution at absolute temperature T given bywhere v0, is the i-component of the velocity of the center of mass of molecule k, andwhere kB is the Boltzmann constant. The temperature T has been chosen to be close to the melting point of the system, by preliminary numerical calculations. In the same way, we assign the magnitude of the angular velocities in eq , whereas the corresponding directions are defined by choosing a random three-dimensional vector.[20] To ensure the mechanical insulation from the exterior, the total momentum of the system has been subtracted. In order to guarantee thermal stability, the velocities are periodically rescaled so that the resulting kinetic energy consistently matches the temperature T within thermal fluctuations. In this way, the system can be considered in contact with a virtual thermal bath. Equation with the initial conditions, eqs –19, has been integrated as described above providing the molecular trajectories

In order to sample a significant region of the phase space, simulations have been performed at different initial conditions and the quantities of interest averaged.

Results and Discussion

We have simulated the trajectories of systems composed of N molecules, ranging from 100 to 256, with many different initial conditions. The dipole and quadrupole have been varied from 0 to 17 D and −15 to 0 B, respectively, thus covering the range of the experimental values of real molecules.[21,22] By preliminary trials, we have succeeded in choosing a set of free force field parameters such that, after a transient initial time t0, molecules aggregate into a single cluster lasting for the whole simulation, for all initial conditions considered.

For each set of trajectories, we have numerically determined the time evolution of the geometric boundary Γ0(t) of the cluster at each time t, t ≥ t0 by reconstruction techniques based on the alpha shapes,[23] and computed the number of atoms of type I (N0I) and II (N0II) on Γ0(t) as well as the number of molecules M0I,II having both atoms on the boundary. These are the key quantities in our study. To analyze the ordering deeper inside the system, we have removed from the data all the molecules having at least one atom on the boundary, evaluated the new boundary Γ1(t), and computed a new set {N1I,N1II, M1I,II}. This procedure is repeated until no more statistically significant internal atoms are present (see Figure ). At the end, a sequence of curves Γ(t), alongside a set of numbers {NI(t),NII(t)),MI,II} l = 0, ...,nmax, is generated at each time step. In the following, we refer to l as the level of boundary, with zero level denoting the outermost curve, also denoted as the interface. The following two-order parameters (OP) have been introduced to quantitatively assess the two possible patterns of orientational order

Figure 1

Snapshot of a cluster comprising 164 molecules (328 atoms) after 7 × 109 MD steps. The system separates into embedded domains delimited by a well-defined set of boundary atoms. As an example, the outermost boundary (l = 0) and an inner boundary (l = 8) are drawn (full line).

Accordingly, O varies from O = 0 [no tangential ordering (TANGO)] to O = 1 (full TANGO) and S varies from S = 0 (no TRANSO) if the same number of type I and type II atoms lie on Γ to S = ±1 where only one type does (full TRANSO in the brackets). For each representation of the system (computer run), the OPs have been averaged over the time t0 < t ≤ tmax, where tmax ranges from 7 × 109 to 1010 MD steps. After a suitable number of runs have been performed, the overall weighed average ⟨O⟩ and ⟨S⟩ with their respective variance have been calculated.

Tangential Ordering

Figure shows the TANGO order parameter ⟨O⟩ versus the radius r of the boundary for a cluster of molecules interacting through van der Waals forces only. It is striking to observe that ⟨O⟩ ≠ 0, ∀l being both dipole and quadrupole moments zero. TANGO appears to be driven by dispersion forces only, originated by intrinsic fluctuations inside the system, which reveals an inherent geometrical nature[24] of TANGO. On increasing the strength of the nonelectrostatic forces, a higher degree of order is observed. This is clearly shown by the upper curve in Figure , which refers to simulations with σ increased by a factor of 1.7 in eq . A more exhaustive analysis would require simulations over a broad range of σ which, even if feasible, is beyond the scope of the present work. It is nevertheless clear that the origin and the magnitude of TANGO is determined by dispersion forces. We remark that the possibility of building superficial ordered structures without electrostatic forces is a peculiarity of TANGO and is excluded in TRANSO, where both multipoles are needed, as will be seen in the next section.

Figure 2

TANGO (both atoms of a molecule on the boundary) resulting from the van der Waals interactions only: the overall order increases for larger ratio σII/σI in eq , that is, for larger ratio of the effective atomic radii. For the same reason, the upper curve is shifted to the right. No electrostatic interaction is explicitly present.

Tangential order is found to be quantitatively independent of the cluster size. Figure also contains the data of a larger cluster (225 molecules). It is noticeable that the data lie on the same curve as the smaller system. This result also rules out that our conclusions might hide some finite size effect.

In the presence of electrostatic forces with one of the two moments being vanished, that is, μ = 0 or Θ = 0, a new important effect is apparent, as shown in Figure : the curve ⟨O0(r)⟩ coincides with ⟨OvdW(r)⟩ of the pure dispersive case, indicating that the dipole–dipole or quadrupole–quadrupole interaction is not able to change the picture. This result further suggests that the dipole and quadrupole moments should enter the functional form of O through their product, in the form

Figure 3

TANGO at zero dipole: all curves for all quadrupole moments overlap, demonstrating the independence of O on the nonvanishing multipole moment. For clarity, only Θ = 0 and Θ = −15 B are plotted. An analogous behavior is found for Θ = 0, μ finite.

A systematic study of the curves O(r) for many values of μ and Θ shows the existence of critical values μc and Θc such that in the rectangle ΓvdW = [0,μc] × [0,|Θc|], all curves ⟨O⟩(r,μ,Θ), plotted as a function of r, overlap, showing TANGO to be insensitive to electrostatics in ΓvdW (data not shown).

As a general rule, TANGO increases over the innermost boundaries, that is, toward regions where the translational invariance tends to be restored. This effect is clearly seen in Figure , where the curves ⟨O(μ)⟩ are plotted at various boundary levels, from the cluster interface l = 0 up to the most internal core which, for a system size N = 100, corresponds to l = 8. The dipole magnitude does not change this general picture, consistently with the geometrical nature of the TANGO. As expected from eq , the same qualitative behavior is found if the quadrupole is varied and the dipole is kept constant (data not shown). Inside the system, the ordering tends to be less affected by the dipole magnitude. This is especially evident on the innermost boundary (upper curve in Figure ), where the order parameter turns out to be basically independent of μ. On the more external layers, the dipole tends to disrupt the order, and, interestingly, this happens above a certain critical value. From eq , it follows that the above results should hold by swapping μ and Θ. In respect to TANGO, we can thus conclude that the effect of the electrostatic interactions is destructive, affecting the outer layers more than the inner ones. At the cluster interface (l = 0), the greatest order is achieved for μΘ = 0. Figure shows this effect in more detail. Deep inside the cluster, the curves ⟨O(r)⟩ at various μ all overlap and, at larger r, separate the pure vdW on the top followed by those with lower μ. The separation occurs at a critical distance of about one-fourth of the total cluster size.

Figure 4

Tangential order parameter vs dipole moment at different levels of boundary, from top to bottom l = 8, ...,0. The lowermost curve (l = 0) corresponds to the cluster interface. Notice the increased ordering inside the cluster at all dipole moments. Quadrupole is kept fixed at Θ = −5 B.

Figure 5

Destructive effect of electrostatics on tangential order affecting the system at all lengths but the inner core.

In the limiting case of strongly interacting molecules, it has been possible to depict a complete phase diagram, for instance, in the (r, μ) plane as shown in Figure . Inside, the cluster is disordered up to a critical distance rc from the center. At larger distances, the destructive effect of the dipole becomes important, and above a critical dipole moment, the system is disordered. Of course, order and disorder are referred to the specific TANGO considered here. Note that, in this context, the phase diagram should be understood bearing in mind that the dipole and quadrupole moments are structural quantities pertinent to the molecules and not environmental parameters as in a conventional one.

Figure 6

Phase diagram of the TANGO in the (r,μ) plane for a system of N = 100 molecules. Above a critical dipole moment, disorder prevails in the outer regions (DIS) in contrast to the internal (W), where van der Waals forces order the system.

Transversal Ordering

In a previous study,[8] it has been established that in a pure LJ system with all atoms having the same value of σ in eq , no TRANSO ordering at any boundary level is achieved. We use this condition in the simulations to “switch off” the LJ interaction, thus disabling any effect on TRANSO, yet keeping both the attractive and repulsive features of the potential which ultimately propel the molecular aggregation. In this way, we can analyze the specific effect of the Coulomb force, excluding any other contribution.

The degree of transversal order at the external–internal interface (l = 0) is shown in Figure as a function of the dipole moment, at various magnitudes of the quadrupole moment. As expected,[8] no ordering is found if either Θ = 0 or μ = 0. In all other cases, the probability of finding a molecule aligned to a direction not tangential to the boundary depends on both the dipole and quadrupole moment. At high values of |Θ|, the functions ⟨S0⟩(μ,Θ) have a very narrow support and approach a delta function. In this regime, a small variation of the dipole moment can change dramatically the order at the interface. At lower values of |Θ|, the transition is less sharp and higher dipole moments are required to achieve comparable values of ⟨S0⟩. The dipole can, however, compensate the lower quadrupole, and for |Θ| as low as about 5 B, complete ordering can still be established, with ⟨S0⟩ = −1.

Figure 7

Order parameter ⟨S0⟩ relative to the outermost boundary (l = 0) as a function of μ, at various quadrupole moments: from top Θ = −0.5, −1, −2, −3, −5, −7.2, −10, and −15 B. Under these settings, ⟨S0⟩ spans from −1 (perfect ordering) to zero (no order). Full lines are fits with Gaussians.

Scaling Behavior at the Interface

From the above considerations, it is clear that the magnitude of the dipole and quadrupole is, to some extent, interchangeable, that is, one can compensate the other as already noticed for TANGO. To examine the point in a more formal way, a model is needed to approximate the functional form of ⟨S0⟩(μ). To this aim, we have fitted the family of Gaussiansto ⟨S⟩(μ,Θ) in the whole range of parameters used in the simulations, with satisfactory results, as shown in Figure . We remark that the choice of a Gaussian as a fitting function is purely empirical with no implication other than its ability to approximate a δ-function. As a result, a single parameter is found to be sufficient to fit the order parameter for all combinations of dipole and quadrupole strengths. This fact again suggests that μ and Θ should not enter S both as independent variables. Going a step ahead, the analysis of the fit parameters shows σ(Θ) to closely follow a power law behaviorwith ν = −1.45 ± 0.02, as shown in Figure . The order parameter can, therefore, be written asthat is, as a function of the sole independent variablewith c being a characteristic length so that λ is a dimensionless quantity. The relevant variable entering the TRANSO order parameter is therefore λ. It sets the scale of the ordering process and ultimately determines the range of validity of our results. This remarkable result is shown by the data collapsing obtained when plotting ⟨S⟩ versus λ reported in Figure . In the limit μ → ∞, Θ finite, and Θ → ∞, μ finite deviations from the scaling are observed. This can be intuitively understood considering such an asymptotic limit as to the lack of one of the multipoles, and we know that such a condition prevents the TRANSO. Such an asymptotic behavior has a merely academic relevance as it can never be reached in a real experiment. Therefore, it will not be considered further in this study. The power law behavior of σ in eq breaks down in the internal layers and should be regarded as a symmetry property of the interface.

Figure 8

Fit parameter σ(Θ) for the data of Figure vs quadrupole moment. The full line represents a power law fitting curve.

Figure 9

Data collapsing of the interface order parameter ⟨S0⟩ on a single master curve after proper rescaling of the multipoles. The error bars on the abscissa are smaller than the symbol size.

Internal Ordering

The order inside the cluster is shown in Figure by the typical set of curves ⟨S⟩(μ) at fixed Θ ≠ 0 and l = 0, ...,8. At all levels l, one can highlight an ordered and a disordered region separated by a transition zone. On increasing μ, the systems are more and more ordered. Such a general behavior is observed at all lengths inside. However, on the internal boundaries, the disordered region becomes wider. Well inside the system, the monotonic behavior of ⟨S(μ)⟩ breaks down, and further increase of the dipole reduces the order due, again, to the predominance of one multipole over the other. By analyzing all curves of this kind in the limit of large Θ, it has been possible to depict the phase diagram in the plane (r,μ) shown in Figure . Deep inside (low r) the cluster is dominated by a disorder and the dipole is not able to change this state irrespective of the strength. This situation changes at a critical distance rc. Above rc, a very unstable region is entered in which a slightly higher dipole generates an ordered system but a further increase would disorder the cluster. This region of instability decreases on increasing r and vanishes at the interface where a critical dipole μc exists such that molecules having μ < μc will generate clusters with a disordered interface and vice versa for μ > μc (see Figure ).

Figure 10

Typical behavior of ⟨S⟩(μ) at fixed Θ = –10 B and at various boundary levels, from left to right, l = 0 (interface) to l = 8 (innermost boundary).

Figure 11

Phase diagram of the transversal order in the (r,μ) plane, at (a) Θ = −15 B in the limit of strongly interacting molecules. At the interface (l = 0), a slight increase of μ is sufficient to generate an ordered system. At lower distances, r order is only possible for values of the dipole moment such that the point (r,μ) lies in the ordered region (ORD). Deep inside the system (l < 6), no ordering is achieved.

From the above considerations, it is clear that the TRANGO and TRANSO are complementary, each may prevail over the other for different regions inside the cluster. In the limit of strongly interacting molecules, on the external interface, TRANSO is the only order possible, whereas TANGO dominates in the interior. This is clearly seen in Figure . As a result, there must be a length where the two orders are of comparable magnitude. In our case, this happens approximately in the center of the cluster and corresponds to the intersection of the two curves shown in Figure .

Figure 12

Comparison of the tangential and TRANSO for strongly interacting particles (Θ = −15 B, μ = 5 D). On the outermost boundary (l = 0), all molecules are oriented with their type I atoms toward the outside (⟨S0⟩ = −1), and no molecule has both atoms on the boundary (⟨O0⟩ = 0). On the innermost boundary, the situation is reversed. Note that at distances r about halfway of the system length, there is the same probability to find a molecule pointing outside and parallel to the corresponding boundary.

Conclusions

We have investigated the formation of ordered structures in molecular clusters by MD simulations. Two basic patterns of order, with molecules oriented tangential or transversal to the boundary, have been considered, which encompass all possible molecular alignment to a boundary, each being measured by an appropriate order parameter. The force field employed has an electrostatic contribution consisting of molecular dipole and quadrupole coupling and a dispersive contribution given by the van der Waals potential. The peculiarity of the force field and the size scale of the observations determine the prevalence of one type of order over the other or their coexistence. For both of them, a phase diagram has been depicted.

From a fundamental standpoint, it is found that TANGO is driven by dispersive forces and is observed even in the absence of electrostatic interactions at all length scales. When present, however, only the dipole–quadrupole interaction is effective and it acts in a destructive manner against TANGO. At the cluster interface, TANGO is greater in the absence of electrostatics and independent of the cluster size. More generally, it is found that the influence of the electrostatic interactions in the ordering process occurs in regions where the translational invariance is broken. In the internal regions where the invariance tends to be restored, electrostatics is irrelevant and the curves ⟨O(r)⟩ overlap for all couples of multipoles considered. For strongly interacting molecules, the general features of TANGO have been summarized in a phase diagram. Inner regions of the cluster are dominated by an order up to a critical length above which systems with a dipole higher than a critical value are disordered. Deep inside the cluster is disordered irrespective of the intensity of the electrostatic forces. Above a critical distance, systems with dipoles higher than a critical value turn out to be ordered.

Transversal order is found to be the result of the broken translational symmetry at the cluster interface and the interaction of multipoles of different parity which, in our second-order simulations, are the molecular dipole and quadrupole. In the absence of electrostatics, the order can still be established, provided the atoms have different σ in the van der Waals potential, so that the required asymmetry can, also in this case, be established. The universality of TRANSO at the interface is revealed by the scaling behavior of the order parameter. This is defined by a scaling exponent and by a characteristic length, which also sets the length scale. In the limit of strongly interacting molecules, a phase diagram can be drawn for TRANSO as well.