Orientation of Electro-optic Chromophores under Poling Conditions: A

J. Phys. Chem. C 2007, 111, 18765-18777

18765

Orientation of Electro-optic Chromophores under Poling Conditions: A Spheroidal Model Harrison L. Rommel and Bruce H. Robinson* Department of Chemistry, UniVersity of Washington, 98195-351700 ReceiVed: May 16, 2007; In Final Form: August 31, 2007

Organic-based chromophores are used as a basic material for electo-optic (EO) devices. The chromophores, to act in an EO material, must be acentrically aligned, which is generally achieved by applying an electric poling field. Calculating, from first principles, the extent of the alignment has proven to be challenging. Therefore, to gain an understanding of the interplay of the fundamental forces that determine alignment, we have chosen to study a coarse-grained model of chromophores: Each chromophore is replaced by a single spheroid, containing a dipole at its center, which interacts with neighboring spheroids only through electrostatic dipolar interactions and a Lennard-Jones potential at the surface of the spheroid. We have found, using NPT and NVT MC simulations that, when the spheroids are allowed to move off-lattice with any orientation, the acentric order from spheroids with varied aspect ratios gives an averaged order which is always lower than that obtained by spheres. The predictions of the order parameters for spheres using off-lattice models are nearly the same as those previously reported for on-lattice models. However, the predictions of the off-lattice models for oblate spheroids (with a 2 to 1 aspect ratio) is in marked contrast to previously described predictions of on-lattice models: The oblate spheroids have approximately half the order of the spheres when off-lattice models are used, as contrasted with having approximately twice the order when using on-lattice models. We find that the acentric order parameters are rather insensitive to the strength of the van der Waals energy over a wide range of energies. The results compare very favorably with measured electro-optic properties of experimental systems.

Introduction Organic chromophores have been developed over the last two decades because they can be used in high-speed electrooptic devices, with large band widths and low drive voltages,1-4 and they are being used as core components of electro-optic and electro-active materials.3,5-11 These chromophores need to be acentrically aligned so that the bulk material has optical activity. Usually, an external electric poling field is applied to drive individual chromophore alignment. The permanent electric dipole moment of each chromophore interacts with the electric poling field and induces a partial alignment of the dipoles, whose average projection onto to the direction of the poling field is 〈cos θ〉. Computing the poling field induced order of chromophores in condensed matter systems has been very challenging.12 Here, we demonstrate that a simple coarse-grained model provides insight into the interactions and is in close agreement with the experimentally determined net order. The experimentally measurable quantity of primary interest is the electro-optic coefficient of the chromophores in bulk. This quantity is directly dependent upon acentric order of the chromophores. The figure of merit is the electro-optic coefficient, r33, of the material. It is related to the phase shift of transmitted light induced by an applied direct current (dc) field in an nonlinear optical (NLO) material (the EO effect). The coefficient r33 is related to the principal components of the molecular hyperpolarizability tensor, β ) β(-ω;0,ω), and the dipole order parameter, 〈cos3 θ〉, by the relation * Corresponding author. E-mail: [email protected].

r33 )

2g(ω)βzzz(-ω;0,ω)FN〈cos3 θ〉 n4

(1)

where βzzz is the projection of the hyperpolarizability tensor onto the dipole axis. This expression is valid when the hyperpolarizability tensor is dominated by a single term that is nearly parallel to the dipole moment, which is often the case.13 Here, FN is the chromophore (or spheroid) number density in units of molecules per cubic centimeter; g(ω) is the product of LorentzOnsager local field factors,8,14 fol and fωl , and the resonance correction for the absorbance of the light. The factors depend on the refractive index of the system and the frequency of the incident light field,

[ ] ( ) ( ( ))

λmax λω g(ω) ) fol (fωl )2 λmax 3 1λω 3-

fol )

(η2 + 2) η2 + 2

fωl )

2

2 2

ηω2 + 2 3

g(ω) is written explicitly in terms of the sample index of refraction, η, dielectric, , maximum absorbance wavelength, λmax, and incident light field8 wavelength, λω, of frequency, ω. For the experimental systems of interest, the quantity g(ω)/η4 is nearly unity. βzzz(-ω;0,ω) is the z,z,z component of the thirdrank molecular hyperpolarizability tensor, responding to light oscillating at frequency ω. The hyperpolarizability is determined

10.1021/jp0738006 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/03/2007

18766 J. Phys. Chem. C, Vol. 111, No. 50, 2007

Rommel and Robinson

by the structural features of the chromophore. The order parameter, 〈cos3 θ〉, is an acentric (odd moment) average tilt angle of the dipole moments of the chromophores from an ordering axis that is co-incident with the external poling field. The cubic dependence of the order parameter arises because of the orientation of each chromophore with respect to (1) the incident light, (2) the applied dc field, and (3) the transmitted light. Another EO coefficient, r13, is important for the EO effect when the geometry of the light field is rotated. The r13 coefficient is related to the average of 〈sin2 θ cos θ〉, so that

r13 〈sin2 θ cos θ〉 〈cos θ〉 - 〈cos3 θ〉 ) ) ) r33 2〈cos3 θ〉 2〈cos3 θ〉

1-

{

2

{

}

〈cos3 θ〉 〈cos θ〉

}

〈cos3 θ〉 〈cos θ〉

The ratio of these two EO coefficients is related to the ratio of the average of the cube and the first power of the cosine. In the limit of noninteracting dipoles and small poling fields, the ratio of these two terms is 3/5. This translates to 1/3 for the ratio of the two EO coefficients. Below, we report 〈cos3 θ〉/ 〈cos θ〉 for the cases of large dipolar interactions. Much work has been done using ellipsoids containing dipoles (Stockmayer fluids) as models for the interaction of molecules (and implicit solvents).15-18 Our interests are to use these models to determine the amount of order that can be developed with an external uniform field, which generates an orienting potential. Such models show characteristics of real systems in that they describe the glassy nature of highly condensed fluids, in which crystalline or ferro-electric behavior are absent. These models provide a method to obtain essential insights into the intermolecular interactions and the net effects of the competition between a poling (ordering) field and the intermolecular electrostatics. Detailed, fully atomistic models of force fields have been used; however, only a few studies have been carried out to understand the acentric order of materials aligned by an external electric (poling) field.12,19,20 Such calculations are numerically intensive, and it is difficult to know when and if equilibrium has been achieved.12,21 Of particular technical difficulty are the many trials and moves needed to achieve equilibrium in such condensed systems. It is of importance to be able to simulate the condensed systems when ellipsoids make up the major part of the molecular system and hence have a large packing density (as a fraction of the packing).22,23 Chromophore densities in typical systems are on the order of 5 to 7 × 1020 molecules/cc, In such dense fluids where the ellipsoid packing fraction, φ, is on the order of 0.20.4, individual moves are difficult to accomplish because of the necessarily small magnitude of allowable trial moves.24-26 Model systems composed of spheroids also serve as a basis for more complex systems that may still be course-grained without the complexity of representing the individual atoms. The use of spheroids as abstractions of molecules is a first step to design molecules of different levels of complexity to find overall structures using detailed force-field interactions. Systems of chromophores on the size of those currently used experimentally27,28 either neat or dispersed in a polymer host cannot be determined with current computational capacity. Coarse-grained models offer the advantage of insight into the effects of simple interactions. Here, we work with a Stockmayerlike fluid: spheroids containing dipoles. The intent is that such a system, with tunable interaction potentials, gives insight into the structure of complex systems and into the gross features

such as order under an external poling field. Moreover, simple models29-32 can be used to test for ferro-electric domains and can be used to determine whether structural features of the underlying molecular systems might induce spontaneous large acentric order in model systems, which may parallel experimental efforts to create highly active electro-optic materials. In a previous paper, we studied the effect of lattice shape on the order parameter for dipoles constrained to rotate on lattice sites.33 In that study, the important finding was that oblate ellipsoids (or discs) were able to find states with order greater than that of spherical ellipsoids. Now, in this paper, we extend our previous simulations to allow the ellipsoids to move offlattice with additional translational degrees of freedom. We find that when this additional translational freedom is added that the greatest order is achieved by spheres: There is no improvement in order under poling for spheroids which have an aspect ratio other than unity. In this work, we report an approach that improves the number of accepted moves without significantly altering the system order, as generated by the poling field and the dipolar interactions. Other methods have been suggested, such as thermal annealing and flipping and rotating methods.29,30 In the present study, the density is kept similar to the experimental systems of around 1 gram/cc; this corresponds to a molecular density around 5 × 1020 to 6 × 1020 molecules/cc. This is most easily done using NVT Monte Carlo (MC) simulations because the density is fixed at the outset (by the size of the spheroid as well as N and V). To allow more rapid equilibration in the NVT MC simulations, the strength of the van der Waals interaction is reduced to where it is no longer strong enough to generate condensation but is large enough to avoid significant spheroidal overlap. To test the effect of this decrease of the van der Waals interaction on the order parameters, we also carried out NPT MC calculations. Under NPT conditions, when varying the van der Waals interactions through the phase transition, from gas to condensed phase, computations demonstrate that the effect on the order parameters is small; these results are consistent with more extensive simulations.15 One conclusion therefore, of this work, is that reduction of the van der Waals interaction does not markedly affect the computed order parameters. Model and Computational Methods Model Potential. An individual chromophore is modeled by an ellipsoid with semi-axes (a,b,c) containing a point dipole at the center of the ellipsoid with a fixed orientation relative to the ellipsoid semi-axes. For the results presented herein, the dipole axis is always aligned with the unique axis (c) of the spheroid. We consider the three cases; spheres (a ) b ) c), oblate spheroids (a ) b > c), and prolate spheroids (a ) b < c). The volume of a spheroid is given by Vspd ) 4π/3abc, The number density (in molecules per cubic centimeter) is FN ) N/V, where V is the volume of the box containing N spheroids. The packing density (fraction) is φ ) NVspd/V ) FN‚Vspd, While oblate off-lattice spheroids can be packed up to a maximum φ ≈ 0.74,26,34 here, we consider packing fractions up to 52%, which is the largest possible fraction for a simple cubic lattice. In comparison, typical chromophore loading is on the order of 25-30% packing density in polymer and up to 45% chromophore in condensed, neat systems.22 The potential energy of the system of spheroids, U, consists of three terms: (1) the interaction of the dipoles of moment, b µ/j , with the poling field, B Epol, (2) the interactions among the dipoles in a dielectric medium, , and (3) the steric interaction of the spheroids with one another via a modified Lennard-Jones

Orientation of Electro-optical Chromophores

J. Phys. Chem. C, Vol. 111, No. 50, 2007 18767

(LJ) interaction potential. These interactions lead to the total system energy, U, as follows:

U ) UPol + UDip + ULJ N

UPol )

N

UDip )

N

ULJ ) 4

{

N′

∑ ∑ i)1 j*i

b µ *i ‚E Bpol ∑ i)1

‚rij3

{ {( ) ( ) ( Roij

N′

∑ ∑ i)1 j*i

i,j

Rij

}

µ/i ‚µ/j - 3(µ/i ‚rˆij)(rˆij‚µ/j )

12

-

Roij

Rij

6

+A

)}}

1 - tanh(f(Rij))

1 f(Rij) ) (|Rij| - 1.12|Roij| - 2) 2

2

first few hundred MC cycles, the overlap of spheroids was shaken out. We also carried out NPT MC calculations to determine the effect of the strength of the LJ interaction parameter on the state of the system and the order parameters in both unpoled15 and poled systems. When computing an ensemble-averaged quantity under constant NPT conditions, the Boltzmann weighting is exp(-1/kBT{U + PV - NkBT ln V}). Here, the energy, U, which is the only energy term in the Boltzmann weighting under NVT simulations, is augmented by the effects of changing the volume of the box, with a constant external pressure, P.36 In the MC algorithm, when the volume is changed by δV, then, the coordinates of the center of each of the spheroids are changed proportionally to reflect the volume change, with no change in the molecular shape or volume. Therefore, U is a function of the volume of the box and is recomputed for the new set of centers of spheroids. The change in the exponential part of the Boltzmann weighting term due to changing the volume by δV is

(V +VδV) ≈ δU + (P - F k T)δV

δU + PδV - NkBT ln The Lennard-Jones energy for each pair of ellipsoids is computed from the ratio (Roij/Rij)2 ) B Roij‚R Boij/R Bij‚R Bij, where B Rij is the vector between the centers of two ellipsoids and B Roij is the vector between the two ellipsoids when they are just touching, when the distance is moved along the vector B Rij. This ratio is the effective contact distance function and is computed from the semi-axes and the relative orientation of the two ellipsoids.35 The LJ energy parameter, i,j, is computed as i,j ) xi‚j, in terms of an energy parameter for each ellipsoid, i.36,37 In the following study, j is taken to be the same value for all spheroids, LJ ) j. In computing the total system energy, the sum over the first index, i, is over all N particles in the box; the second index, j, is summed over all other spheroids within a minimum cutoff distance, subject to the use of re-entrant boundary conditions and the minimum image convention.33,36,38 Monte Carlo. We carried out NVT MC simulations on the three different spheroid systems at a variety of densities. The MC moves are performed on the spheroids sequentially. The system consists of 512 spheroids per box. The individual spheroids are permitted to rotate and to translate (for off-lattice calculations only). Since we are interested in the bulk response in the absence of any external shape effects, all runs were performed with a cubic simulation cell with re-entrant conditions. Systems of 512 spheroids and larger were sufficient to eliminate any finite size effects because of the periodic boundary.33 For each MC sweep through the lattice, each spheroid is rotated and translated by a random increment (as a rigid body transformation), and the move is accepted or rejected according to the Boltzmann weighting and the Metropolis criterion.36 The size of the rotation (and translation) is determined during a training session so that approximately half of the moves are accepted.36 The magnitude of the translations and rotations typically scale with the system density. The order parameters are compiled from the final 1000 accepted MC moves for 100 independent MC trajectories of length 5000 moves at each volume (or density). Systems typically required on the order of 100-200 MC steps to achieve stability, as determined by the fluctuations about the mean energy and order parameter. As a test to determine whether the moves were still active, the field was removed, and the systems returned to a centrosymmetric order within 200 MC steps. All of the offlattice calculations begin with the spheroids placed on the lattice sites of a simple cubic lattice with random orientation, allowing for possible interpenetration of neighboring spheroids. In the

N B

The only difference between the NPT and the NVT MC computations is that after all possible translations and rotations of each of the individual spheroids are completed (as done in an NVT MC) a change in the volume is made and the energy is recomputed and tested against the Metropolis acceptance criterion. In both NVT and NPT calculations, the order parameters were measured as a function of the strength of the LJ energy parameter, LJ, which was varied from 0.1kBT to 2kBT. In NPT calculations, it was found that the system condensed at higher LJ, consistent with previous work. This gives the domain of LJ over which the system order remains constant. Results In most of the following studies, the poling field is fixed at Epol ) 150 V/µm; the dipole moment is µ* ) 24 D; the dielectric is ) 4, and the LJ energy parameter LJ ) 0.12kT. The temperature is T ) 350 K, corresponding to kT ≈ 0.05 picoergs. The magnitude of the dimensionless poling parameter is f ) µ*Epol/kT ) 2.4, The aspect ratio of the prolate and oblate spheroids (of a:c) are 1:2 and 2:1, respectively. In all cases, the volume of the spheroids remains constant regardless of the aspect ratio, and the radius of the sphere is a ) b ) c ) 5.57 Å, The volume of a spheroid then is 0.0724 × 10-20cc, The dipole moment and molecular volume correspond to those needed to simulate ONLO molecules (such as YLD156 which is similar to FTC38,39) in a dendritic environment. The poling field is typical of the experimental maximum; the temperature is near the experimental poling temperature for the dendritic material containing ONLO chromophores, and the dielectric constant is near that determined experimentally.27 Figure 1 shows a schematic of a box containing either prolate (left) or oblate (right) spheroids. The field is oriented vertically; each spheroid contains a dipole (represented by the black arrow) at the center oriented co-incident with the unique axis of the spheroid. The angle between the unique axis of the spheroid and the poling field is θi. The order parameters are then N cosn θi, where n ) 1, 2, 3, or computed as 〈cosn θ〉 ) 1/N∑i)1 4 and summed over the N molecules in the box. The odd values of n ) 1 or 3 are the two acentric order parameters of interest. Figure 2 shows the loading parameter for the three types of spheroids (the prolate, spherical, and oblate) in an on-lattice


Figure 1. Orientation of 1:2 prolate spheroids (left) and 2:1 oblate ellipsoids (right) and the angle (θi) with respect to the poling field; the poling field is co-incident with the Z axis.

calculation. The lattice is a simple cubic lattice distorted to match the aspect ratio of the spheroids. This is the same deformation used in other studies to demonstrate the effect of spheroid distortion on order parameters.40-42 The effect of lattice choice (among face centered cubic (FCC), body centered cubic

Rommel and Robinson (BCC), and hexagonal close packed (HCP)) makes small difference to the resultant order parameter, as long as the deformation by an affine transformation of the lattice matches the spheroid deformation, as previously presented.33 For comparison, the low-density limiting values for one- and threedimensional systems are shown as straight solid lines. For the one-dimensional case in which the dipoles may point only either with or opposite to the poling field, at low density, the order parameter for both n ) 1 and n ) 3 is 〈cosn θ〉 ) tan h(f) ∼ f, and for the three-dimensional case, 〈cosn θ〉 ) Ln(f) ∼ f/(2n + 1), where Ln(f) is the nth Langevin function43 and the approximate value is obtained only in the low f limit. In this onlattice model, the oblate ellipsoids achieve an order that exceeds that of the low-density limit for the three-dimensional system and is about twofold larger than that obtained for spheres (under large density conditions). Figure 3 (top) shows the first acentric order parameter (n ) 1), 〈cos θ〉 for spheres, on-lattice and off-lattice as a function of the density of the spheres, FN. The order parameter for different lattice types (labeled on the figure) is also shown, demonstrating that the nature of the lattice does not lead to dramatically different results for spheres. Moreover, the order parameter for spheres, when allowed to move off-lattice, while somewhat larger than any lattice model, is nearly the same as for the on-lattice spheres. As expected, for the dilute limiting case regardless of model, the same value of the order parameter is obtained, and it is the same as the Langevin result for dilute dipoles. Similar results are shown for the 〈cos3 θ〉 order parameter, showing that the effect of moving off-lattice does

Figure 2. Loading parameter, LP ≡ FN〈cos3 θ〉, as a function of density for 1:2 prolate ellipsoids (diamonds), spheres (circles), and 2:1 oblate ellipsoids (squares). The spheroids are constrained to remain on the lattice. The lattice points (corresponding to an affine deformation of a simple cubic lattice) are co-incident with the ellipsoid aspect ratios. Solid lines denote the theoretical predictions of the Ising limit, 〈cos3 θ〉 ) tan h (f), and the dilute Langevin limit, 〈cos3 θ〉 ) f/5.33,38,43 The dipole moment is µ* ) 26 D, the poling field Epol ) 150 V/µm, the dielectric ) 4, and the temperature is 373 K. Therefore, the dimensionless poling parameter is f ) µ*E/kBT ) 2.5.



Figure 3. Top: Comparison of the order parameter, 〈cosn θ〉, n ) 1, as a function of the spheroid density FN, for the simple cubic (blue squares), body-centered (green, no marker), and face-centered cubic lattice (red diamonds) with the off-lattice model (black circles) for dipolar spheres of radius r ) 5.5 Å. Also shown is a comparison of the n ) 3 order parameter for on-lattice and off-lattice models both using spheres. The dipole moment is µ* ) 24 D, the poling field is Epol ) 150 V/µm, the dielectric is ) 4, the temperature is 350 K, and f = 2.5. The upper limit on the abscissa, FN ) 7.25 × 1020 molecules/cc, corresponds to that of closest packing for spheres on a simple cubic lattice, φ ) 0.52. The error bars are the standard error for the parent distribution from 105 MC trajectories compiled from 100 independent simulations. Bottom: The centrosymmetric order parameter, 〈P2〉 ) 1/2(3〈cos2 θ〉 - 1), as a function of number density computed from the same trajectories as used to compute the data in the top figure.

not significantly change either the n ) 1 or n ) 3 acentric order parameter for spheres and that off-lattice spheres obtain only slightly larger order at high density. Over the entire range of spheroid density, the ratio of the two acentric order parameters is nearly constant, bounded by

L3(f) L1(f)

) 0.65 g

〈cos3 θ〉 g0.6 〈cos θ〉

At low density, FN < 2 × 1020, the ratio tends to the Langevin function ratio, and at high density, the ratio tends to 0.6. The value of 0.6 is the ratio of the Langevin functions in the perturbation limit of f , 1. Figure 3 (bottom) shows the centrosymmetric order parameter

1 〈P2(cos θ)〉 ) (3〈cos2 θ〉 - 1) 2


Rommel and Robinson

Figure 4. Order parameter, 〈cos3 θ〉, dependence on the strength of the poling field, Epol, at low, FN ) 0.05 × 1020/cc, and high, FN ) 7.25 × 1020/cc, number density, for spheres (blue circles) and prolate spheroids (red diamonds) computed by NVT-MC. The analytic, Langevin result, L3(f) (dashed line), is overlaid and is close to the low-density result. Straight lines are the best-fit lines proportional to the poling field for the two types of spheres. The first-order parameter, 〈cos θ〉, (not shown) is (5/3)〈cos3 θ〉.

as a function of the number density for the three different spheroid types. In the absence of a poling field, 〈P2(cos θ)〉 ) 0. 〈P2(cos θ)〉 is rather insensitive to the aspect ratio of the spheroids. A comparison of the top and bottom curves shows that the centrosymmetric order parameter, 〈P2(cos θ)〉, is about 1.6 fold smaller than 〈cos3 θ〉 at low density and 2- to 4-fold smaller at high density. At low density, the ratio of 〈P2(cos θ)〉 to 〈cos3 θ〉 is due to the ratio of Langevin functions, L2(f) and L3(f). At high density, the poling-induced order is a small perturbation to the other energies driving the system to no net order. As such, at high density, 〈P2(cos θ)〉 is proportional to the poling parameter squared, f2, whereas the acentric order parameters are proportional to the poling parameter to the first power, f. Therefore, centrosymmetric order parameters may retain a smaller poling-induced order at high-density relative to their asymmetric counterparts.44 Figure 4 shows the loading parameter, LP ) FN‚〈cos3 θ〉, dependence on the strength of the poling field, Epol, at both low FN ) 0.05 × 1020/cc and high FN ) 5.0 × 1020/cc number density for the off-lattice spheroids. The functional form of the loading parameter at low density follows that of the Langevin function, L3(f). At high-density, however, the order parameter and hence the loading parameter are linear in the poling field strength. At high density, the poling field energy is small relative to the dipolar and van der Waals energy terms and can be considered a perturbation. The first-order expansion of this perturbation term shows that the acentric order should be proportional to the poling field. These results hold for all three spheroid types.33,45 Figure 5 shows the loading parameter, LP ) FN‚〈cos3 θ〉, for the three off-lattice spheroids as a function of the spheroid number density, FN. This figure may be directly compared with Figure 2, in which the calculations are nearly identical here except the spheroids are allowed to move off lattice. The ratio of 〈cos3 θ〉/〈cos θ〉 (not shown here) is the same as for the on-

lattice models (Figure 2), that is ∼0.6. This result applies to all three shaped spheroids. The error bars on each of the calculations represent the standard error from the set of 500 distinct MC equilibrium runs at each different density. Because the error bars are the standard error for the parent distribution, the errors in the mean values are negligible. The solid lines are only to aid the eye in following the trends. At high number density, FN > 4 × 1020, the loading parameter for the oblate spheroids is about half that of spheres, and the loading parameter for the prolate spheroids is intermediate between the spheres and the oblate spheroids. This result is the central difference between on-lattice and off-lattice models. As illustrated in Figure 6, the neighbor pairwise correlation function, g(r), illustrates the influence of spheroid shape on order. g(r) is a continuous function of the distance to a neighbor, averaged over all centers. Because we use re-entrant conditions (minimum image convention), the maximum distance that can be covered is half the length of the simulation cell. Therefore, r is restricted 0 er e L/2, where V ) L3. The distances are divided on a uniform grid into units of ∆r, so that the continuous function may be replaced by the discrete histogram gm ) g(rm ) m∆r). The r dimension is divided so that rm ) m∆r and the integer index, m, runs from 1 to M, where M∆r ) (1/2)L. Nl(rm) is the number of spheroids over the distance interval (rm - ∆r) er erm from any center. The average number of molecules, averaged over C centers is

N h (rm) )

1

C

∑Nl(rm)

C l)1

Associated with the average number of spheroids from the center is the pairwise correlation function defined as



Figure 5. The loading parameter, LP ) FN‚〈cos3 θ〉, for off-lattice model simulations, as a function of spheroid number density, FN, for 2:1 oblate ellipsoids (black triangles, c ) 3.5 Å), 1:2 prolate ellipsoids (red diamonds, c ) 8.84 Å), and dipolar spheres (blue circles, c ) 5.57 Å). Other parameters and error bars are computed using the same prescription as Figure 3. All three systems approach the Langevin limit at low density.

Figure 6. Pair wise correlation function for a high-density case, FN ) 7.25 × 1020 molecules/cc, for the spheres (blue circles), oblate (black triangles), and prolate (red diamonds) spheroids. The minimum distances of closest approach are 7, 8.8, and 11 Å for the oblate, prolate, and spherical cases, respectively. Also shown are the integrated number densities (inset), which are the average number of neighbors out to the distance specified on the abscissa.


Rommel and Robinson

Figure 7. Histogram of the interior angle, γ, between each spheroid and its nearest neighbors for a low density (FN ) 0.5 × 1020 molecules/cc), top, and a high density (FN ) 7.25 × 1020 molecules/cc), bottom, case for spheres (blue circles), oblate (black triangle) and prolate spheroids (red diamonds). The strong peak near 〈cos γ〉 ∼1 indicates head-to-tail ordering within the system. All other parameters are the same as Figure 5. The solid line in the low-density plot is the independent sphere, normalized Boltzmann distribution for the relative angle.46

g(rm) )

N h (rm) FN4πrm ∆r 2

Therefore, the average total number of neighboring molecules p out to distance r ) rm must be NT(rp) ) ∑m)1 N h (rm) ) p 2 ∑m)1g(rm)FN4πrm ∆r. As r becomes large, the molecular density should approach the average density:

NT(r) rf∞ 4π 3 r 3

FN ) lim

The first “peaks” occur at distances comparable to the minimum (hard sphere) distances of 7, 9, and 11 Å for the oblate, prolate, and spherical spheroids, respectively. The spheroid overlap results from the competition between the repulsive part of the Lennard-Jones energy and the dipolar interaction energy (discussed later). The close approach for the oblate case, which occurs at all densities, indicates the preferred head-to-tail, pairwise alignment of these particular spheroids.

At low density for the oblate spheroid case (data not shown), there is no second peak in the pairwise distribution function, g(r), indicating that there is no clustering of next-nearest neighbors. At high density (FN ) 7.25 × 1020/cc), g(r) indicates dimer formation at an optimal distance of ∼6 Å for the oblate spheroids and further generalized clustering of about 5 neighbors as evidenced by the second (smaller peak) around 12 Å. The total number of neighbors (on average) within 12 Å is around 5 for all three shaped spheroids (see inset). Figure 7 shows the histogram of spheroids as a function of the relative orientation angle, γ, of the neighbors, from the same MC data depicted in Figure 5. The low-density distributions are all very much like the expected low-density Boltzmann distribution (shown as a solid line).46 The number of spheroids near the relative angle of cos γ ∼ 1 indicates a preference for head-to-tail orientations for the oblate spheroids, and a preference against head-to-tail orientations for the prolate spheroids, relative to the simple spheres. At high density (FN ) 7.25 × 1020/cc), there is minimal correlation of the relative angle, because the overall order with respect to the poling field is so much lower (see Figure 3). One can average over as many



Figure 8. “Snapshot” of the oblate spheroids at low (top) and high (bottom) density from NVT calculations. The densities are the same as those of Figure 7. The color-coding indicates the direction of the dipole oriented in the unique direction of the spheroid.

nearest neighbors as desired with no significant change in the cos γ distribution function. There is no strong preference of orientation of nearest neighbors at any density except for the well-defined (head-to-tail) orientation of oblate spheroids, in the region of cos γ ∼ 1. Figure 8 shows a snapshot (in the equilibrated form) of the oblate ellipsoids at the same densities for which the histograms (Figure 7) are shown from NVT calculations. The snapshot shows the head-to-tail local ordering observed in the histograms in Figure 7. The high-density example has a packing density of φ ) 0.51, which is the same that one has for a simple cubic lattice. As can be seen in the low-density snapshot, the oblate spheroids tend to arrange head to tail in a partially condensed coexistence with the gas phase. The preferential orientation and pairing carry over to the more condensed phase (as illustrated in the histogram of Figure 7). This then leads to lower acentric order parameters for the oblate spheroid case. Figure 9 shows the results using NPT MC simulations of spheres with (triangles) and without (circles) dipoles. In these simulations, the external pressure is held constant (at 2 atm),

and the volume of the box is allowed to change. The simulations were performed in a manner similar to the NVT cases; the system was allowed to come to equilibrium with 2500 accepted MC moves followed by 2500 equilibrium accepted moves averaged to obtain the system density. The error bars represent the standard distribution of densities obtained by 50 independent MC trajectories. At a low LJ energy parameter, the density is consistent with an ideal gas. When LJ > kBT, the density asymptotes to FN ≈ 5 × 1020 molecules/cc or φ ≈ 0.35 for the liquid phase with no dipolar interactions and asymptotes to slightly less when the dipoles are included. This density (and related packing fraction), when LJ > kBT, indicates a condensed (liquid) phase. The computations shown in the figure were done with no poling field. The inclusion of the dipolar interactions (triangles) causes the system to condense at a somewhat lower LJ energy parameter and obtain a slightly reduced density in the condensed phase. In general, the inclusion of a poling field (not shown) does not markedly alter the shape of the curve. The acentric order parameters, determined from NPT MC computations at LJ ) 1.0 × kBT and FN ) 4.5 ×


Rommel and Robinson

Figure 9. Equilibrium system density determined from the NPT MC method,36 at an external pressure of 2 atm, as a function of the LJ energy parameter (LJ), with dipoles (triangles) and without dipoles (circles) present, and no poling field. Otherwise, the same parameters are used in these computations as in the previous NVT computations (Figures 3 and 5).

1020/cc are consistent with results obtained from NVT computations at the same density when LJ ) 0.1 × kBT. Figure 10 shows the average energies per molecule as a function of the strength of the Lennard-Jones interaction, LJ, for both the dipole-dipole energies (triangles), Eµµ, and the Lennard-Jones energies (circles), ELJ, in the NPT MC simulations of spheres. The error bars represent the standard error of the parent distribution. The results show that the dipole-dipole energy remains essentially constant in the range of the condensed phase, LJ g 0.6 × kBT. Condensation to the liquid phase appears to be relatively independent of electrostatic energy; condensation is driven primarily by the attractive part of the LJ interactions. Even at the lowest feasible LJ (in the gas phase), the dipolar interactions are relatively constant and only decrease minimally due to increased spheroid overlap, increasing favorable pairwise dipolar interactions. By contrast, the LJ energy monotonically decreases as the LJ energy parameter increases, even after the spheroids have formed a condensed phase. Discussion Off-lattice models allow for types of packing that permit a more realistic description of the order of a system of chromophores under poling conditions than on-lattice models. The poling strengths, dipolar energies, and the molecular dimensions are realistic for the experimental chromophore systems being modeled in these studies. However, to ensure that enough freedom is given to the system of spheroids, the LJ energy parameter has been reduced to one that is about 1 order of magnitude smaller than values known for similar, realistic systems.47 Only when the LJ energy is kept low (less than kBT) is it possible to make MC moves in NVT simulations that reach

equilibrium in a reasonable number of steps, always less than 500. The LJ energy, when using a low LJ energy parameter, is a continuous function of intermolecular distance that effectively retains a semihard-sphere picture of the spheroids. The strength of the attractive part of the LJ energy is not a significant part of this model. The LJ energy is small enough to allow spheroidal movement but large enough to keep spheroids from significant overlap and unrealistically large dipolar interactions. In studying more complicated systems where multiply connected ellipsoids are used to model dendrimeric structures,27 it was found that low values for the LJ energy parameter were required to prevent trapping, which would result in low acceptance of MC moves. This was one of the primary motivations to study the effect of varying the LJ energy in simple systems. As the nearest neighbor correlation function shows, the low-energy parameter does still generate a realistic system, with only a small amount of spheroid interpenetration. This is the motivation to explore the utility of using small values of the LJ energy parameters. One obtains the same order parameter (within error) from NVT MC calculations using a small LJ energy parameter as one obtains from NPT MC simulations with a large LJ energy parameter for identical densities. Here, we now consider the verification and validation of the computational methods used to arrive at this conclusion. First, under dilute concentrations the MC methods (both NVT and NPT) arrive at the same results as analytic theory. Using NPT methods and a low LJ energy parameter, LJ, the density is the same as that expected for an ideal gas. The phase transition observed with increasing LJ is consistent with the phase transition described by Smit and Leeuwen15 and Weis and Levesque48 using a similar model. Our results are consistent with their finding that only the



Figure 10. Average energies per molecule as a function of the LJ energy parameter (LJ) for an NPT MC simulation of dipolar spheres. Both the dipole-dipole energy (triangles) and the Lennard-Jones energy (circles) are shown. Parameters are the same as previous figures.

attractive part of the LJ energy leads to a condensed system. The dipolar interactions by themselves do not appear to lead to a condensed system. Second, a variety of tests were performed to ensure that the training of step sizes gave an approximately 50% acceptance ratio in most runs. In general, the magnitude of rotational and translational trial moves decreases with increasing density. Checking for the presence of sufficiently large fluctuations in the order parameters during the course of MC trajectories, on the order of the error bars depicted in Figures 5 and 8, ensured that the systems were adequately sampling the equilibrium phase space and not trapped in metastable energy minima. Individual MC runs with different random number seeds all gave similar fluctuations and equilibrium order parameters. Third, the internal order parameter and the neighbor distance parameters indicate that at high density one obtains multiple nearest neighbors, indicative of condensed phases. However, at all densities, the oblate spheroids showed a tendency to have an additional single nearest neighbor at roughly twice the minor semi-axis distance, indicative of dimer formation, as can be seen in Figure 8. Fourth, in the models developed, the strengths of the poling and dipolar fields were on the order expected from experimental systems.22,23,28,49 The order parameters obtained in these studies can be used to compute the EO coefficient. For example, the order parameter is 〈cos3 θ〉 ≈ 0.15 at the density FN ) 6 × 1020/cc, see Figure 5.50 The experiments done on a trimer of YLD156 using dendritic synthesis yielded a response of r33/Epol ) 1.42 ( 0.04 (nm/V)2 for fields up to a poling field of Epol e 150 V/µ.27 The hyperpolarizability, βzzz ) 2800 × 10-30esu, is estimated from HRS experiments and is consistent with first principles computation.13 The simulation of the loading parameter using the off-lattice prolate-spheroid model and eq 1, g(ω)/η4 ≈ 1, yields a predicted r33theory ) 210 pm/V using a poling field of Epol ) 150 V/µm. The EO

coefficient, r33, is proportional to the poling field strength. Therefore, r33theory/Epol ) 210/150 ) 1.4 (nm/V)2, in agreement with the experimental results. A variety of initial conditions were tried, all of which gave the same order parameters once equilibrium had been achieved. For example, the final equilibrated order parameters were independent of initial orientation: random, ferroelectric, and antiferroelectric initial conditions all converged to the same average order parameter. The shape of the box containing 512 spheroids was changed (the aspect ratio was varied up to a factor of 2): The order parameters were insensitive to these shape changes. Moreover, the poling field was turned on and off to test that the system was never locked into any particular state. Whenever the field was removed, the acentric order parameters, as well as 〈P2(cos θ)〉, returned to zero (results not shown). The combination of information obtained from the pair wise correlation function and the relative orientation (data in Figures 6 and 7) combine to suggest a grouping of 5 or 6 neighbors for the prolate and spherical cases, which were centrosymmetric in the absence of the poling field. The oblate spheroids show dimerization, going head to tail, as well as the clustering observed for the prolate and spherical cases. These results are consistent with the previous work on dipolar spheres which showed that the LJ dispersion term is necessary to go from a vapor to a liquid state.15 Our NPT calculations are carried out through the gas-liquid phase transition predicted by van Leeuwen and Smit.15 A realistic LJ energy, sufficient to condense the spheroids, causes the MC algorithm to slow down in convergence. This is a well-known problem with MC calculations of condensed phases.24,36 However, by employing a gas-phase LJ parameter (i.e., LJ ∼ 0.1kT), the NVT MC simulations produce accurate order parameters within a reason-

18776 J. Phys. Chem. C, Vol. 111, No. 50, 2007 able computation time, avoiding the problem of computational slowing. Conclusions The dramatic change in poling-induced order due to shape for the on-lattice model entirely disappears when the spheroids are allowed to move off lattice. The basic result is that there is a small, twofold, dependence on overall shape of the spheroids (or chromophores). More importantly, off-lattice models show that there is no way to improve on the order of spherical objects in the absence of lattice organizing forces. As an example, the oblate spheroids, according to the off-lattice model, at high number density, have an order that is twofold worse than that of spheres. The parameters used in this study are typical of the chromophores now being developed to make materials with large, that is, r33 > 300 pm/V, electro-optic coefficients. The present theoretical results compare exceedingly well with experimental systems. In experimental studies, YLD156, as part of a dendrimer system, gave an EO response27 of r33/Epol ) 1.4 (nm/V)2, which is in excellent agreement with our present result. This same system has also been studied for its centrosymmetric order, in which the measured order51 was found to be 〈P2〉/Epol ) 8 × 10-4 µm/V. This compares well to the present theoretical result, Figure 3B, of 〈P2〉theory/Epol ) 5 × 10-4 µm/V. The off-lattice models provide three additional degrees of freedom for each molecule. Therefore, we take the off-lattice computations to be superior to any on-lattice computational model in that an artificial constraint has been relaxed. The important results are threefold. First, the results for spheres do not change between on- and off-lattice models. This is important because it illustrates the validity and consistency of the methods. Packing is quite similar among different lattices and so offlattice is not expected to be markedly different for spheres for off- or on-lattice models. However for spheroids, the effect of lattice geometry does cause marked differences in order parameters depending on the assumptions of the lattice structure. Therefore, for nonspherical spheroids, it is not surprising that the off-lattice models give markedly different results. The second point is that the prolate spheroids consistently show lower order for both on- and off-lattice models. However, when the lattice constraint is removed, the ordering of the prolate spheroids is not far below that of the spheres, being reduced by around 7080%. Third, the oblate spheroids are particularly important because the on-lattice calculations predicted that one could obtain an order better than that achievable in the low-density limit (of simple Langevin theory). This would be truly remarkable if possible. However, the off-lattice calculations show that the order expected from oblate spheroids is even worse than that achievable from prolate spheroids and is about 50% that of spheres. This relative reduction in order is due to the presence of dimer formation as evidenced in the data of Figure 6. While on-lattice models of oblate spheroids predict the formation of strong ferroelectric chains that enhance net ordering, in contrast, off-lattice simulations show that head-to-tail ordering occurs only at the dimer level. Beyond dimer pairs, other geometries compete with an extended head-to-tail ordering. A dimer of two oblate spheroids has the approximate shape of a sphere (albeit with a larger dipole), and spheres do not readily form dimers. This behavior is a consequence of the additional translational freedom of the chromophores and suggests that spheres have the proper shape to achieve optimal order. Oblate spheroids can achieve high order only by invoking additional organizational forces above and beyond those generated by a poling field and intermolecular dipolar interactions.

Rommel and Robinson Acknowledgment. The authors acknowledge the support of the STC-MDITR Program of the National Science Foundation (DMR0120967 and NSF-DMR-0092380). Support from the Air Force Office of Scientific Research under AFOSR-(F4962003-1-0110-P000), as well as from the DARPA MORPH program Phase I ((N) 14-04-10094), is also gratefully acknowledged. Any opinions, findings, conclusions, or recommendations expressed in this work are those of the authors and do not necessarily reflect the views of these funding agencies. We thank Dr. Robert D. Nielsen for many thoughtful comments on this manuscript. References and Notes (1) Lee, M.; Katz, H. E.; Erben, C.; Gill, D. M.; Gopalan, P.; Heber, J. D.; McGee, D. J. Science 2002, 298, 1401. (2) Marder, S. R.; Cheng, L.-T.; Tiemann, B. G.; Friedli, A. C.; Blanchard-Desce, M.; Perry, J. W.; Skindhoj, J. Science 1994, 263, 511. (3) Shi, Y.; Zhang, C.; Zhang, H.; Bechtel, J. H.; Dalton, L. R.; Robinson, B.; Steier, W. H. Science 2000, 288, 119. (4) Luo, J.; Liu, S.; Haller, M.; Liu, L.; Ma, H.; Jen, A. K.-Y. AdV. Mater. 2002, 14, 1763. (5) Dalton, L. D. AdV. Polym. Sci. 2002, 158, 1. (6) Dalton, L. R.; Harper, A. W.; Ghosn, R.; Steier, W. H.; Ziari, M.; Fetterman, H.; Shi, Y.; Mustacich, R. V.; Jen, A. K. Y.; Shea, K. J. Chem. Mater. 1995, 7, 1060. (7) Zyss, J. Molecular Nonlinear Optics; Academic Press: San Diego, 1994. (8) Prasad, P. N.; Williams, D. J. Introduction to Nonlinear Optical Effects in Molecules and Polymers; John Wiley and Sons: New York, 1991. (9) Dalton, L. R. Pure Appl. Chem. 2004, 76, 1421. (10) Jen, A. K.-Y.; Neilsen, R.; Robinson, B.; Steier, W. H.; Dalton, L. Mater. Res. Soc. Symp. Proc. 2002, 708, 153. (11) Kim, T.; Kang, J.; Luo, J.; Jang, S.; Ka, J.; Tucker, N.; Benedict, J.; Dalton, L. R.; Gray, T.; Overney, R.; Park, D.; Herman, W.; Jen, A. K. Y. J. Am. Chem. Soc. 2007, 129, 488. (12) Leahy-Hoppa, M. R.; Cunningham, P. D.; French, J. A.; Hayden, L. M. J. Phys. Chem. A 2006, 110, 5792. (13) Isborn, C. M.; Leclercq, A.; Vila, F. D.; Dalton, L. R.; Bre’das, J. L.; Eichinger, B. E.; Robinson, B. H. J. Phys. Chem. A 2007, 111, 1319. (14) Onsager, L. J. Am. Chem. Soc. 1936, 58, 1486. (15) van Leeuwen, M. E.; Smit, B. Phys. ReV. Lett. 1993, 71, 3991. (16) Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1957, 27, 1208. (17) Wertheim, M. S. J. Chem. Phys. 1971, 55, 4921. (18) Hansen, J.; McDonald, I. Theory of Simple Liquids, 2nd ed.; Academic Press: London, 1986. (19) Kim, W. K.; Hayden, L. M. J. Chem. Phys. 1999, 111, 5212. (20) Makowsa-Jankusik, M.; Reis, H.; Papdopoulos, M. G.; Economou, I. G.; Zacharopoulos, N. J. Phys. Chem. B 2004, 2004, 588. (21) Tu, Y.; Luo, Y.; Agren, H. J. Phys. Chem. B 2006, 110, 8971. (22) Sullivan, P. A. Chapter 2: Theory Guided Design and Molecular Engineering of Organic Materials for Enhanced 2nd-Order Nonlinear Optical Properties. Ph.D. Thesis, University of Washington, 2006. (23) Sullivan, P. A.; Akelaitis, A. J. P.; Lee, S. K.; McGrew, G.; Lee, S. K.; Choi, D. H.; Dalton, L. R. Chem. Mater. 2006, 18, 344. (24) Landau, D. P.; Binder, K. A Guide to Monte Carlo Simulations in Statistical Physics; Cambridge University Press: Cambridge, 2005. (25) Ferrenberg, A. M.; Landau, D. P.; Binder, K. J. Stat. Phys. 1990, 63, 867. (26) Donev, A.; Cisse, I.; Sachs, D.; Variano, E. A.; Stillinger, F. H.; Connelly, R.; Torquato, S.; Chaikin, P. M. Science 2004, 303, 990. (27) Sullivan, P.; Rommel, H.; Liao, Y.; Olbricht, B. C.; Akelaitis, A. J. P.; Firestone, K. A.; Kang, J.-W.; Luo, J.; Davies, J. A.; Choi, D. H.; Eichinger, B. E.; Reid, P. J.; Chen, A.; Jen, A. K.-Y.; Robinson, B. H.; Dalton, L. R. J. Am. Chem. Soc. 2007, 129, 7523. (28) Liao, Y.; Anderson, C. A.; Sullivan, P. A.; Akelaitis, A. J. P.; Robinson, B. H.; Dalton, L. R. Chem. Mater. 2006, 18, 1062. (29) Ayton, G.; Patey, G. Phys. ReV. Lett. 1996, 76, 239. (30) Ayton, G.; Wei, D. Q.; Patey, G. N. Phys. ReV. E 1997, 55, 447. (31) Weis, J. J. J. Chem. Phys. 2005, 123, 044503. (32) Camp, P. J.; Patey, G. Phys. ReV. E 1999, 60, 4280. (33) Nielsen, R. D.; Rommel, H. L.; Robinson, B. H. J. Phys. Chem. B 2004, 108, 8659. (34) Sloane, N. J. A. Nature 1998, 395, 435. (35) Perram, J. W.; Wertheim, M. S. J. Comp. Phys. 1985, 58, 409. (36) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (37) Hildebrand, J. H.; Scott, R. L. Regular Solutions; Prentice Hall: Engle Cliffs, 1962.

Orientation of Electro-optical Chromophores (38) Robinson, B. H.; Dalton, L. R. J. Phys. Chem. A 2000, 104, 4785. (39) Liu, S.; Haller, M. A.; Hong, M.; Dalton, L. R.; Jang, S.-H.; Jen, A. K.-Y. AdV. Mater. 2003, 15, 603. (40) Pereverzev, Y. V.; Prezhdo, O. V. Phys. ReV. E 2000, 62, 8324. (41) Pereverzev, Y. V.; Prezhdo, O. V.; Dalton, L. R. J. Chem. Phys. 2002, 117, 3354. (42) Pereverzev, Y. V.; Prezhdo, O. V.; Dalton, L. R. Chem. Phys. Lett. 2001, 340, 328. (43) Piekara, A. Proc. R. Soc.: London A 1939, 172, 360. (44) Wallace, P. M.; Sluss, D. R.; Dalton, L. R.; Robinson, B. H.; Reid, P. J. J. Phys. Chem. B 2006, 110, 75. (45) Dalton, L. R.; Harper, A. W.; Robinson, B. H. Proc. Natl. Acac. Sci. U.S.A. 1997, 94, 4842.

J. Phys. Chem. C, Vol. 111, No. 50, 2007 18777 (46) Nielsen, R. D. The Relative Angle Distribution Function in the Langevin Theory of Dilute Dipoles. ArXiv.org: http://arxiv.org/list/ physics.class-ph/recent, (accessed March 5, 2007). (47) Evans, D. J.; Watts, R. O. Mol. Phys. 1976, 32, 93. (48) Weis, J. J.; Levesque, D. J. Chem. Phys. 2006, 125, 34504. (49) Dalton, L.; Robinson, B. H.; Jen, A.-K.; Ried, P.; Eichinger, B.; Sullivan, P. A.; Akelaitis, A. J. P.; Bale, D.; Haller, M.; Luo, J.; Liu, S.; Liao, Y.; Firestone, K. A.; Bhatambrekar, N.; Bhattacharjee, S.; Sinness, J.; Hammond, S.; Buker, N.; Snoeberger, R.; Lingwood, M.; Rommel, H.; Amend, J.; Jang, S.-H.; Chen, A.; Steier, W. H. Proc. SPIE-Int Soc. Opt. Eng. 2005, 5935, 593502/1. (50) Kim, T.-D.; Kang, J.-W.; Luo, J.; Jang, S.-H.; Ka, J.-W.; Tucker, N.; Benedict, J. B.; Dalton, L. R.; Gray, T.; Overney, R. M.; Park, D. H.; Herman, W. N.; Jen, A. K.-Y. J. Am. Chem. Soc. 2007, 129, 488. (51) Olbricht, B. C.; Sullivan, P. A.; Davies, J. A.; Takayesu, J.; Reid, P. J.; Dalton, L. R.; Robinson, B. H. J. Appl. Phys. 2007, submitted.

Orientation of Electro-optic Chromophores under Poling Conditions: A

Recommend Documents