How Ethylene Glycol Chains Enhance the Dielectric Constant of Organic Semiconductors: Molecular Origin and Frequency Dependence.

Incorporating ethylene glycols (EGs) into organic semiconductors has become the prominent strategy to increase their dielectric constant. However, EG's contribution to the dielectric constant is due to nuclear relaxations, and therefore, its relevance for various organic electronic applications depends on the time scale of these relaxations, which remains unknown. In this work, by means of a new computational protocol based on polarizable molecular dynamics simulations, the time- and frequency-dependent dielectric constant of a representative fullerene derivative with EG side chains is predicted, the origin of its unusually high dielectric constant is explained, and design suggestions are made to further increase it. Finally, a dielectric relaxation time of ∼1 ns is extracted which suggests that EGs may be too slow to reduce the Coulombic screening in organic photovoltaics but are definitely fast enough for organic thermoelectrics with much lower charge carrier velocities.

Introduction

The excitonic nature of organic semiconductors (OSCs), namely, the generation of bound electron–hole pairs instead of free charges upon photon absorption, has been attributed to their low dielectric constant.[1−4] Koster et al.[5] argued that high dielectric constant OSCs can reduce the exciton binding energy and consequently increase the efficiency of organic photovoltaics (OPVs), provided that the energy offset required to enable charge transfer between acceptor and donor is minimized. Since then, addition of ethylene glycol (EG) side chains to fullerene derivatives,[6−8] small molecules,[9,10] and polymers[11−13] has become the prominent strategy for enhancing the dielectric constant of OSCs.[4] However, increased dielectric constants have not resulted in higher power conversion efficiencies so far,[4,14] and the effects on the exciton binding energy have not been reported. On the other hand, for organic thermoelectrics (OTEs), the use of EGs has shown to improve the thermal stability, doping efficiency, and power factors.[15−19] This effect is attributed to the polar environment of EGs,[19] likely related to the increased dielectric constant.

An important aspect that often receives little to no attention is that the dielectric constant enhancement due to EGs has already completely vanished at the high-frequency limit (∼1 THz),[9,11] meaning that the enhancement is not electronic but due to nuclear relaxations. It should be noted that this is different for nonexcitonic silicon solar cells where the full dielectric constant of ∼12 is electronic.[20] This raises the question whether the slower EGs can provide additional Coulombic screening in OSC applications, as their response (i.e., their reorientation) would need to be faster than the mobility of the charge carriers. In other words, if the charge carriers change their environment faster than the environment can respond, no additional screening can be expected from the nuclear relaxation of the environment.

It has been shown, using the time domain reflectometry technique, that for pure liquid EGs, the transition frequency for the nuclear dielectric contribution is 1 GHz[21] and that this contribution completely vanishes in the solid state, which leaves EGs with a dielectric constant below 3.5.[22] Therefore, it is of interest to understand if and how this contribution persists in solid state organic electronics and identify its time scale. Since impedance spectroscopy has been the only technique used in previous experimental work so far for organic electronics with EGs[4,6−13] and because it can only probe up to the megahertz regime,[23,24] where the EG contribution is still fully active, the transition frequency for the dielectric contribution of EGs, i.e., the time scale of their response, is yet unknown to date.

In the present work, we study the dielectric constant of a fulleropyrrolidine with a single EG side chain named PTEG-1 (Figure ) which has been recently synthesized and shown to have a static dielectric constant of 5.7 ± 0.2.[6] As this is higher than the solid state dielectric constant of both C60 (∼4)[25] and EGs (∼3.5),[22] there is clear evidence for a synergistic effect that is not fully understood yet. We employ our new computational protocol (see Methods section) to calculate both its time- and frequency- dependent dielectric constant. Having an atomistic resolution, we are able to pinpoint the different dielectric contributions (electronic, dipolar, induced) to different fragments of the molecule and identify the molecular response that is causing the unusually high dielectric constant. Our calculations yield accurate static and electronic dielectric constants and allow us to predict the dielectric transition frequency and the time scale of the nuclear relaxations which are then used to investigate the relevance of the dielectric constant increase due to EGs for OSC applications.

Figure 1

PTEG-1 molecular structure and a sample snapshot from a simulation box. Different molecular fragments, as used in this work, are highlighted in different colors.

Results and Discussion

In Figure , the dielectric constant versus both time (blue) and frequency (red) is presented. Static and electronic dielectric constants of 5.97 ± 0.03 and 3.41 are calculated, respectively, in good agreement with both the experimental[6] static (5.7 ± 0.2) and the theoretical[25] (periodic coupled perturbed density functional theory) electronic (3.3) dielectric constants. The dielectric relaxation time of the system is computed to be 0.95 ± 0.05 ns, which corresponds to a transition frequency of 0.17 ± 0.01 GHz, which we predict to be the frequency at which this relaxation can be observed experimentally. This is approximately five times slower than the pure liquid transition frequency of triethylene glycol,[21] likely due to the reduced flexibility in the solid phase. Dependence of this transition frequency on different molecular features and different morphologies—potentially resulting from different processing conditions[26]—is currently under investigation. We anticipate that our computational protocol can be highly beneficial to help the endeavor of engineering molecules with faster dielectric responses.

Figure 2

Computed dielectric constant versus time (blue) and frequency (red) at 25 °C. ϵ0 and ϵ refer to the static and electronic dielectric constants, for which the experimental[6] and theoretical[25] references are 5.7 ± 0.2 and 3.3, respectively. Imaginary part of the dielectric constant is shown in Figure S3.

In order for OSCs to benefit from the static dielectric constant, the charge carriers would need to change their environment at a slower rate than the dielectric relaxation time of 0.95 ns. Using the electron mobility and internal electric field in OPV devices one could approximate the electron hopping rate as ∼2 ns–1.[27] It has also been argued that the actual charge carrier motion in OPV devices is orders of magnitude faster than what would be expected based on their mobilities[28−30] due to the nonequilibrium nature of OPVs, which would result in correspondingly higher electron-hopping rates. This suggests that the nuclear response of EGs could be too slow to influence the Coulombic screening for OPVs. On the other hand, for OTEs, the much smaller electric field in the devices results in an approximate electron-hopping rate of 2 μs–1,[31] which is orders of magnitude slower than the dielectric relaxation time. Moreover, in such doped organic semiconductors the ionized dopants are essentially static, allowing EGs to permanently screen their charges. These results are in line with the current performance of EG containing OSCs in these respective fields: While no improvement to OPV efficiencies has been made due to the inclusion of EG side chains,[4,14] for OTEs, inclusion of EGs has been shown to enhance device performance.[15−19] This indicates that it is crucial to look more carefully at the transition frequency and that the static dielectric constant should not be taken by default as the effective dielectric constant for OSC applications.

The computational protocol presented here further allows, as shown in Figure , a clear-cut decomposition of the dielectric constant (see Methods) into contributions of molecular fragments (as defined in Figure ) and dielectric processes (electronic, dipolar, induced). The electronic contribution is due to the response of the electrons before any nuclear movement occurs. Then the nuclear response is split into dipolar and induced contributions: the former is the dipolar alignment of the molecule due to the partial charges on each atom, while the latter is the additional electronic polarization due to the dipolar alignments, i.e., a contribution that is coupled between the electronic and the nuclear parts. The importance of a polarizable molecular dynamics (MD) simulation becomes apparent at this point, as from a classical “fixed-charge” MD simulation[32] only the dipolar contribution could have been obtained. The results show that the electronic contribution is dominated by C60 as it has a highly polarizable electron cloud. This shows that increasing the size of the side chain would result in an overall reduced electronic dielectric constant, as was also concluded in previous work.[25] The dipolar contribution can be almost fully attributed to the EG groups, implying a significant reorientation of the EG dipoles, which is clearly shown later in Figure . The induced contribution is shared almost equally between the C60 and the EG fragments, suggesting a strong dielectric enhancement due to favorable interactions between the highly polarizable C60 and the highly dipolar EG, which we argue as an important reason for the synergistic increase of PTEG-1’s dielectric constant (5.7) compared to its C60 (∼4) and EG (∼3.5) components in the solid state. The total contributions show that the group connecting C60 to EG, named “connection”, provides very little contribution overall, meaning that minimizing its size can be a design rule for increased dielectric constants. While similar total contributions for C60 and EG are observed, considering that EG is about one-half the size of C60, its contribution per volume is much more significant. However, since different fragments are responsible for maximizing the electronic and nuclear contributions, the trade-off and the need for focusing on the more relevant contribution for OSCs becomes even more apparent.

Figure 3

Decomposition of the dielectric constant into contributions from molecular fragments and from dielectric processes. Total of each bar is shown in bold type. Difference of 1 between the contributions shown here and the static and electronic dielectric constants from Figure corresponds to the vacuum dielectric constant (see eq ).

Figure 4

P1 order parameters for the individual ethylene glycol (COC) groups. COC vector is defined as r⃗C – r⃗C (see also inset).

We now further analyze whether the dipolar contribution coming from the EGs is indeed due their significant reorientation in response to the electric field. To this end, the order parameter P1 (cos θ) is calculated for each of the EGs as a function of time, where θ is the angle between the direction of the applied field and the COC vector, as shown in Figure . An order parameter P1 = 0 means random orientation of the EGs, and P1 = 1 means perfect alignment with the direction of the field. Before the electric field is applied (t = 0), all EGs are randomly oriented as one would expect in an amorphous system. With the sudden application of the field, EGs orient in the direction of the field at a similar time scale as the dielectric relaxation time (Figure ), which, together with the results from Figure , shows that the alignment of the EG groups is indeed responsible for the high dielectric constant. Moreover, it can be seen that the EG group directly connected to the benzene ring aligns the least with the field, while that at the end of the side chain aligns the most. This suggests a dependence of the EG contribution either on the total length of the chain, i.e., longer chains give higher contributions toward their end, or simply on the position of the EG within the chain, i.e., the terminal EG has increased flexibility. We suggest that in the former case EG chains branching toward the end, and for the latter case multiple short chains could result in increased dielectric contributions.

Finally, to identify what particular feature of EGs allows them to easily align with the field even in the condensed phase, we show the energy profile (Figure ) for their two distinct torsions within the EG fragments: OC–CO (blue) and CO–CC (red), which we obtain by performing an inverse Boltzmann analysis (see Methods) to their distributions throughout the simulations. Then the transition rates between the minima are approximated by transition state theory. For both torsions, a region of approximately 270° with three minima can be seen where the torsional barrier is always lower than 8 kJ/mol. Having such a region allows the reorientation of EG at a picosecond time scale, which clearly is the reason for the high flexibility of EGs.

Figure 5

Energy profile of the OC–CO (blue) and CO–CC (red) torsions averaged over all such torsions within the EG fragments, all of the molecules, and all of the simulations. Numbers accompanied by the arrows indicate the transition rate over a barrier in the specified direction.

Conclusions

In summary, we outlined a computational protocol that can predict the time- and frequency-dependent dielectric constant of organic solids with high accuracy. We demonstrated this using the PTEG-1 molecule, which contains EG side chains that are known to significantly enhance the static dielectric constant even in the solid state. We showed that this enhancement occurs by the alignment of EGs with the electric field, which is in turn made possible by their low torsional barrier. We made several design suggestions to maximize the static dielectric constant, such as minimizing the size of the group connecting C60 to EGs and using shorter or branched EG chains. Moreover, we identified the transition frequency and the corresponding dielectric response time of PTEG-1 as ∼0.2 GHz and ∼1 ns, respectively, which was yet unknown to date. Due to their very different charge carrier velocities, we argued that this response is fast enough to be fully benefitted by OTEs, while it may be too slow to provide additional Coulombic screening in OPVs, which is also in agreement with the current performance of EGs in these fields. Therefore, the static dielectric constant, often obtained by impedance spectroscopy, is not necessarily the effective dielectric constant for OPVs, and more efforts should be made on decreasing the response time of these nuclear contributions or on increasing the electronic dielectric constant instead of the static one. We believe that the computational protocol described in this work can be highly beneficial for these efforts.

Methods

Computation of the Time- and Frequency-Dependent Dielectric Constant

While the most common and convenient approach to compute the dielectric constant from MD simulations is by monitoring the fluctuations of the dipole moment in an equilibrium simulation,[33−36] this approach requires unreasonably long simulation times in the case of solid systems, making it unsuitable for such applications. Moreover, the fluctuation method is not able to capture the electronic contribution to the dielectric constant even when used in conjunction with polarizable force fields as no such fluctuations occur during the simulation.[37,38]

The external field method,[33,34,38,39] as used in this work, allows for much shorter simulation times, and even though it requires a higher number of simulations, these are embarrassingly parallel. It also makes it possible to obtain the electronic dielectric constant when used in conjunction with polarizable force fields. Moreover, having a nonequilibrium simulation with an applied field allows us to directly look at the molecular response to the electric field, similar to how molecules would respond to the field generated by charge carriers in OSCs. In this method, an electric field Eext is applied in the direction i and the dipole moment μ(t) is monitored as a function of time. Then the time-dependent dielectric constant ϵ(t) is given bywhere V is the volume of the simulation box and μinit is the initial dipole moment before the electric field is applied. The applied field method also has its own challenges: The strength of the applied field must be chosen carefully for each system since too small values make it difficult to distinguish the dipolar response from the statistical noise and at too large values the linear relationship between the dipole moment and the field strength no longer holds.[38,40] Furthermore, the usage of this method is much less straightforward than the dipole fluctuation method, requiring various scripts and steps, which makes it less accessible for the nonexpert user.

For this work, a computational protocol that provides an easy pipeline to compute the dielectric constant of any system with the applied field method was developed, provided that the user has an appropriate force field for the system of interest. The applied field methodology used in Riniker et al.[38] was combined with a fitting procedure to extrapolate the simulation to longer time scales, and by performing a Fourier transform to this fit, the frequency-dependent dielectric constant was obtained. The dielectric response, ϵ(t), as shown in Figure , has a sharp initial increase followed by a single negative exponential decay. The exponential decay part was then fit to such a functionwhere λ is the difference in ϵ(t) between the start of the fit and the converged value and τ is the dielectric relaxation time. Having λ as a fitting parameter allowed capturing the last few percent of the dielectric response without needing excessive simulation times due to the slow converging nature of the exponential decay function as also seen in Figure . A similar approach was also followed by Van der Spoel et al.[41] for the dipole fluctuation method where they extrapolated the dipole autocorrelation function to convergence. The authors provide the necessary scripts for the use of this protocol with the GROMACS[42] software in the Supporting Information.

The step-by-step procedure is as follows: (1) Multiple condensed phase morphologies are generated by performing high-pressure “squeezing” simulations (see Supporting Information for details); (2) for each of the morphologies, an equilibrium simulation is run and a snapshot is taken at given intervals; (3) for each of the snapshots three new simulations are started with an applied electric field in the x, y, or z direction; (4) the time-dependent dielectric constant ϵ(t) is computed using eq ; (5) the average ϵ(t) from all of the simulations is extrapolated to longer time scales using eq ; (6) this fit is then Fourier transformed to obtain the frequency-dependent dielectric constant; (7) the dielectric constant is decomposed into dielectric processes and molecular fragments of interest, as further explained in the next section.

In this work, we used a Drude-based polarizable force field[43] that is based on our newly developed Q-Force procedure in which force field parameters are derived from quantum mechanical calculations in an automated way (see Supporting Information). Three amorphous PTEG-1 morphologies were generated, and for each of them, 46 snapshots are taken with 100 ps intervals after an initial relaxation of 1 ns, resulting in 414 simulations over which all of the results are averaged over.

Decomposition of the Dielectric Constant into Molecular Fragments and Dielectric Processes

The electronic dielectric constant was obtained at t = 0 of the simulation by looking at the difference between the dipole moment before and after the application of the electric field which is caused by the displacement of the Drude particles. The dipolar contribution was obtained by recalculating the dipole moment with the Drude particles put back on top of their corresponding atom. This removed all of the contribution due to Drude particles, resulting in the contribution that is solely due the reorientation of dipoles which are originated due to partial charges on the atoms. The induced contribution was then what remains after subtracting the electronic and dipolar contributions from the total dielectric contribution, which can be visualized as further polarization of Drude particles as a function of time, i.e., the coupling between the electronic and the dipolar contributions.

Decomposition of the dielectric contribution into molecular fragments was done by considering the dipole moments of the fragments of interest instead of the whole molecule. It is important to note that since the individual fragments are not necessarily uncharged, their dipole moment becomes origin dependent. However, for calculation of the dielectric constant, the quantity of interest is the derivative of the dipole moment with respect to the applied electric field, which again becomes origin independent. The MDAnalysis library[44,45] was used to apply the transformations mentioned above.

Torsional Free Energy Barriers and Transition Rates

Free energy profiles for the CO−CC and OC−CO torsions were obtained through an inverse Boltzmann analysis of the torsional distributions that were averaged over all simulations and all torsions of that type. The relative energy of the torsion with angle α (Eα) is given bywhere R is the gas constant, T is the temperature, nα is the number of occurrences of the angle α, and N is the total number of data points. Transition rates were approximated using transition state theory assuming a two-state free energy difference usingwhere kb, h, and R are the Boltzmann, Planck, and gas constants, respectively, T is the temperature, and ΔG is the free energy difference between the minimum and the corresponding transition state. It is important to note that this is a crude approximation, neglecting the coupling between the torsions, and is aimed at only giving an approximate time scale.

Statistics and the Error Margin

The statistical error was calculated from the standard error of the mean (SEM) usingwhere σ is the standard deviation and n is the number of simulations. As seen in Figure S4, the standard deviation is 0.47 with the dielectric constant varying between 5 and 7. This means that it is very important to have sufficient simulations to obtain meaningful results. With the 414 simulations that we have, the statistical error margin becomes ±0.02. Then we approximated the error due to using a single-exponential fit by varying the starting time of the fit (varied between 100 and 600 ps), for which we determined an error of ±0.01. We also obtained the error margin for the transition frequency and the relaxation time with this method. For the presented dielectric constant, we summed both of these error margins, which resulted in a total error margin of ±0.03. It is important to note that there are additional sources of error, such as the choice of force field, but it is nontrivial to quantify these errors.

Molecular Dynamics Run Parameters

The double-precision version of GROMACS 2018.x[42] software was used for all simulations. The equations of motion were integrated using a leapfrog algorithm with a time step of 2 fs. The cutoff for Lennard–Jones interactions was 1.2 nm. NPT ensemble was used in all simulations: The temperature was set to 298 K and the pressure to 1 bar, unless stated otherwise. The Berendsen[46] thermostat (coupling parameter = 1 ps) and barostat (coupling parameter = 5 ps, compressibility = 4.5 × 10–5 bar–1) was used. The electrostatic interactions were treated by the particle mesh Ewald (PME) method.[47] The strength of the applied electric field was 0.25711 V/nm (0.0005 au). Coordinates along the trajectories were written up to the fifth decimal due to the very small displacement of Drude particles. Further details can be found in the MD parameter file in the Supporting Information.