Non-negligible Outer-Shell Reorganization Energy for Charge Transfer in Nonpolar Systems

Many charge-transporting molecular systems function as ordered or disordered arrays of solid state materials composed of nonpolar (or weakly polar) molecules. Due to low dielectric constants for nonpolar systems, it is common to ignore the effects of outer-shell reorganization energy (λout). However, ignoring λout has not been properly supported and it can severely impact predictions and insights derived. Here, we estimate λout by two means: from experimental ultraviolet photoelectron spectra, in which vibronic progression in these spectra can be fitted with the widths of peaks determining the low-frequency component in reorganization energy, regarded to be closely associated with λout, and from molecular dynamic (MD) simulation of nonpolar molecules, in which disorder or fluctuation statistics for energies of charged molecules are calculated. An upper bound for λout was obtained as 505 and 549 meV for crystalline anthracene (140 K) and pentacene (50 K), respectively, by fitting of experimental data, and 212 and 170 meV, respectively, from MD simulations. These values are comparable to the inner-sphere reorganization energy (λin) arising from intramolecular vibration. With corresponding spectral density functions calculated, we found that λout is influenced both by low- and high-frequency dynamics, in which the former arises from constrained translational and rotational motions of surrounding molecules. In an amorphous state, about half of the λout’s obtained are from high-frequency components, which is quite different from the conventional polar solvation. Moreover, crystalline systems exhibit super-Ohmic spectral density, whereas amorphous systems are sub-Ohmic.


INTRODUCTION
Electron transfer (ET) theory 1 has been an important foundation for understanding charge transfer processes in many systems.One important application area is in chargetransporting materials, such as organic light-emitting diodes, solar cells, and field-effect transistors, etc. 2−8 With their advantages of low cost, low mass, and flexibility, organic molecular semiconductors have attracted much attention. 9In the Marcus theory, ET rate can expressed as where V is the electronic coupling between initial and final states, and ΔE 0 is the zero-point energy difference of the two states.The reorganization energy λ is a critical factor in the Marcus theory, as it is tightly linked to the activation energy needed for electron transfer.−17 λ comprises two components: the inner-sphere λ in , arising from intramolecular degrees of freedom, and the outer-shell, λ out , from changes in the surrounding environment.The former can be estimated from standard quantum chemistry calculations 18,19 in which the energy difference of two reactant and product structures in the same state is calculated, with the electronic state being that for either reactant, or product, in the calculation.The latter, λ out , can be calculated using a simple expression 1,20 = + e r r R ( ) 1 2 2 opt s i k j j j j j y { z z z z z i k j j j j j j y { z z z z z z (2)   where the two fragments involved in the ET, the donor D and the acceptor A are modeled as spherical cavities with radii r D and r A , with their center-to-center distance R, embedded in a dielectric medium characterized by static and optical dielectric constants ε s and ε opt , respectively.According to eq 2, nonpolar solvents have very similar ε s and ε opt values, leading to vanishing λ out , e.g.20−30 meV for benzothieno-benzothiophene (BTBT) with long alkyl side-chains, 21 and 33−42 meV for fullerene. 22−26 Values of λ out can be inferred from experiments, and most of them are on the order of 100 meV, similar to many λ in values reported.In ref 27., an experimental work, λ out 310 and 140 meV were fitted from data obtained with electron attachment.λ out can also be estimated from dynamic Stokes shifts in a number of nonpolar solvents, which can reach 150 meV.Interestingly, a correlation for the Stokes shift value to effective quadrupole moments was observed. 28By fitting emission spectra with a semiclassical Marcus equation, the reorganization energy derived from the width of the vibronic progression is largely contributed by λ out .In a porphyrin-fullerene system in benzene, it was estimated as 130 meV., 29 and 80−532 meV for various donors in fullerene-based organic solar cell systems. 30heoretical estimates of λ out for nonpolar solvents are seen with various approaches.Based on dipolar solvation theories with structural factors considered, 31,32 the program SolvMol yielded solvent reorganizations of 2.58 and 0.24 meV for metaand para-phenyl-acetylene-bridged carbazole-naphthalimide in weakly polar toluene. 33With the density structural factor of solvents considered, an expression for the induced component in λ in response to structural changes, λ ind was developed, 34 and it is estimated as 100−300 meV for nonpolar solvents, a range that fits experimental observations.Even though the dipolar contribution does not exist in nonpolar solvents, a work focusing on quadrupole solvation was developed, 35 and values of 165−205 meV were reported for the reorganization energy of benzene, a range that is similar to the magnitude of commonly reported λ in for organic molecules exhibiting good semiconducting potential (In ref 23., λ in values reported are 138 and 102 meV, for tetracene and pentacene, respectively.)Therefore, with dielectric polarization theory, both the density structural factor (structural density fluctuation) in electric induction, and the quadrupole moments of solvent can lead to significant values in λ out In the Marcus theory, the reaction coordinate is a collective quantity over many degrees of freedom.Free energy curves ΔG of the initial and final states are assumed to be parabolic, where λ is the curvature.Therefore, the width of a thermal distribution gives rise to the reorganization energy.With this approach, values of 10−100 meV were reported for nonpolar systems. 36,37−46 This is usually rationalized with the low dielectric constants in nonpolar environments in those studies, together with the expression in eq 2. Works employing molecular modeling, such as hybrid quantum mechanical/ molecular mechanical (QM/MM) schemes and polarizable force fields, λ out were calculated with optimization in several common organic molecular crystals, but only 0.98−7.81meV were obtained. 47,48In these works, the energy of one of the charged states was calculated at different optimized structures, a typical method for determining λ in .Doing so for λ out by including the surrounding molecules may sound like a reasonable generalization, but statistical contributions were ignored.
While direct comparisons of computational characterizations of charge transporting behaviors with experimental measurements is not entirely feasible yet, we believe that including contributions of λ out in theoretical models is important.Therefore, careful examination of the magnitude of λ out in typical nonpolar systems is necessary.
In the present work, we estimate λ out from experimental data, as well as direct simulation for nonpolar molecules.We studied ethylene, anthracene, and pentacene in their solid states, with ethylene being a small system suitable for benchmarking.−51 The experimental bandwidth of ultraviolet photoelectron spectra (UPS) was employed to estimate λ out .MD simulation was used to account for fluctuation in the interaction energy of a cationic molecule embedded in the system, from which λ out was calculated.Both results show values larger than 100 meV, values that cannot easily be ignored.We also studied spectral density functions and their breakdown in frequencies, fitting for the Ohmic model.We further discuss the role of λ out in ET rates.

THEORY AND METHOD
2.1.Extracting Reorganization Energy from Ultraviolet Photoelectron Spectra.In UPS, an incident highenergy photon ionizes a molecule, ejecting an electron.The energy of the cationic state of the system can be inferred from the kinetic energy of the electron ejected.In this way, UPS spectra show transitions from neutral molecules to cations, with possible vibronic progressions.Since it is half of a hole transfer process, these spectra contain information about the reorganization energy of hole transfer in the medium employed in the experiments.
2.1.1.Theoretical Ground.Based on Marcus-Levich-Dogonadze-Jortner (or Jortner's) theory, for a system with a high-frequency mode, the electron transfer rate can be expressed as 52−54 where V is the electronic coupling between initial and final states, and ΔE 0 is the zero-point energy difference of the two states.In this widely recognized formulation, it has been assumed that the reorganization energy can be split into lowand high-frequency contributions, using thermal energy k B T as the dividing line.λ s is the low-frequency contribution in the reorganization energy, and for the high-frequency contribution, it is usually written with just one representative (effective) mode, with ω v being the frequency of the mode, and S, the Huang−Rhys (HR) factor, with the high-frequency part of reorganization energy λ v as, The original multimode expression can be found at refs 52, 55.
Reducing multiple high-frequency modes to an effective mode representing the vibronic progression effect has also been discussed. 52,56Multimode with quantum statistical treatment for ET rate has been an active area of theoretical works, and in general, systems without a clear vibronic progression in their density of states can also be treated. 57Molecules with strong vibronic progression spectra usually have more than one active, high-frequency mode.As long as they are with similar frequencies, an effective frequency and an effective HR factor can be defined, 56,58 and eq 4 with such effective factors gives the reorganization energy of these high-frequency modes.
With vibronic progression in a transition between ionic and neutral states, as in eq 3, the HR factor S for a high-frequency vibration that has the largest coupling to the ionization can be inferred. 599][50][51]60 By taking the standard deviation σ of the progressive band, together with Gaussian terms in eq 3, we have In practice one can infer the full width at half-maximum (fwhm) easily from experimental spectra, and use the following relationship for a Gaussian function, 9][50][51]60 Therefore, in the present work, we took the fwhm from experimental literature and used it in the theoretical model composed of Gaussian bands for estimating λ s eqs 5 and 6 as illustrated in Figure 1.

Partitioning in λ and
Possible Discrepancies between λ out and λ s .Illustrative depiction for two different ways to partition the reorganization energy λ are included in Figure 1a.Traditional inner-and outer-shell separation has λ out arising from intermolecular interactions, comprised of solute− solvent, or system-bath interaction, and λ in mainly arising from intramolecular vibrational degrees of freedom.On the other hand, there exist active high-frequency modes that contribute to the progression of spectra, and the remaining low-frequency modes, and their corresponding portion of reorganization energy are modeled as λ s and λ v in fitting the photoelectron spectra.Such division allows a simplified classical treatment for the low-frequency mode and the quantum effect of highfrequency modes can be included, leading to Jortner's rate expression (eq 3) or other multimode variations. 52,56,57t is desirable to assume that the two partitionings of λ are highly correlated, i.e., it is likely that the intermolecular interaction is dominated by low-frequency motions, and the reorganization energy from intramolecular vibrational modes is mainly high-frequency, and out s in v (7)   If so, estimating λ s from UPS spectra yields λ out directly.−63 In other words, there might be a high-frequency contribution in λ out that belongs to λ v , and such a contribution is quite small for the polar solvents reported.On the other hand, low-frequency modes, such as ring-distortion, as depicted in Figure 1a, contribute to λ in due to its intramolecular nature, but it is possible to contribute to λ s with its low frequency.Therefore, rigorously speaking, we have two terms in the discrepancy between λ out and λ s : Another possible discrepancy of this approach is that widths of these progression bands can be broadened by various frequencies of high-frequency modes.The λ s from fitting the width of the vibronic bands might be overestimated.
In the present work, we estimated λ s 's from experimental UPS spectra and cautiously compare them to calculated λ out .In doing so we simply assumed that these two additional terms in eq 8 are small or they could largely be canceled.Therefore, values of λ s from fitting UPS spectra in the literature can be treated as an upper-bound for λ out .We note that it is not a fully rigorous approach as further details are remained to be explored.

Estimating λ out with Molecular Dynamics Simulation. 2.2.1. Energetic Disorder and Its Relationship to λ out .
Reorganization energy is the energy required in changing the external coordinates from those in equilibrium with the initial state, to those of the final state, without changing the electronic configuration.For the outer part, λ out , it is the energy difference arising from system-environment interaction.For example, if we have a hole-transfer event in two identical molecules where we keep the ordering of the first and the second molecules M. The reorganization energy is, where E(A|B) denotes the (free) energy of state A in the system equilibrated under state B. Since outer-shell reorganization energy involves changes in the environment, in this work, we consider the long-range, electrostatic interaction between the system and surrounding environment.The interaction energy of the additional charge on one molecule and the environment ΔE int (t) is defined as, where E int neutral (t) is the interaction energy of a target molecule and all other molecules, with the target molecule in its neutral state, and E int charged (t) being the same interaction energy but with the target molecule in its charged state.After extracting snapshots from MD simulation, we calculate the energy difference of a (arbitrarily chosen) target molecule of its charged and neutral states.In other words, ΔE int (t) is the MD realization of intermolecular, electrostatic contribution of E(M + |M) − E(M|M), whose statistics lead to λ out .
In principle, twice the ensemble average for ΔE int is the outer-shell reorganization energy value we are seeking, λ out .Since we have MD trajectories, the dynamic aspect of λ out can be inferred as well.Equivalently, the reorganization energy can be evaluated using the fluctuation−dissipation theorem. 61,64he autocorrelation function (ACF) C(t) can be defined as where δX(t) ≡ΔE int (t) − ⟨ΔE int ⟩ and ⟨ΔE int ⟩ is the mean value of ΔE int (t) in each trajectory.Therefore, with Fourier transform of ACF { } C t ( ) the power spectrum I(ω) can be obtained as, In the present work, we employ the spectral density J(ω) to describe the fluctuation of electrostatic energies.Under the classical limit, it can be derived from I(ω): where we have included a factor of 2 since we were only calculating half of the ET process (M → M + ) in the simulated ΔE int (t).λ out is calculated as, 61,64 To evaluate the quality of theoretical results of λ out from MD trajectories, the standard deviation (STD) of the mean value of the autocorrelation function C(t) is employed.With eqs 14 to (17), and by accounting the factor of 2 in eq 12, λ out , can also be expressed as, In the MD simulation, we include 20,000 ethylene, or 2,000 anthracene or pentacene molecules in the simulated box.Each trajectory is recorded at 1 fs time steps, for 18 ps.With each MD trajectory (time-evolution of a box of equilibrated, neutral molecules), we collect more than 300 trajectories of ΔE int (t), by placing a positive charge on a different molecule.Molecules selected to carry positive charges were separated, as we required at least 10 Å for ethylene, or 7 Å for antracene and pentacene, between the molecules selected.
To properly take ensemble averages, we have repeated the accumulation of data from multiple independent trajectories of MD simulation, with 300 sampled ΔE int (t) from each trajectory.The uncertainty of λ out values was calculated from the standard deviation (SD) of averaged C(0) in eq 18, until the standard deviation of the averaged λ out is 1% or less of their mean values.In so doing, we accumulated more than 3,700 MD trajectories in each case and the final number of trajectories varied case by case.Under this setting, the averaged long-time correlation function values (C(t) at 9 ps) are of the order 10 −4 eV 2 .
2.2.2.Other Details of Molecular Dynamic Simulation.Force field parameters were generated with Q-Chem and Q-Force packages. 65,66In quantum chemistry calculations, ground state geometry optimization and vibrational analyses were performed with long-range corrected density functional LC-BLYP/DZ* with dispersion correction 67,68 included.Longrange correction parameters, μ, employed for ethylene, anthracene, and pentacene are 0.41, 0.2, and 0.2 Bohr −1 , respectively.Q-Force was employed to generate force-field parameters by fitting with force generated by the ab initio method.
In turning a neutral molecule to a cation, we calculate charges from the electrostatic potential on a grid (CHELPG), 69 and add charge differences of its neutral ground state and its cationic state to the force field of the cation.In this calculation, we took the globally optimized structure calculated at LC-BLYP/DZ* with dispersion correction.In order to capture the system-bath interaction for λ out , we did not update the intramolecular force field parameters, such as equilibrium bond lengths and bond angles when a molecule is in the cationic state, as they contribute to λ in .
We employed OpenMM to simulate molecular dynamics. 70o reach an equilibrated state, the system was prepared with an Andersen thermostat and a Monte Carlo barostat, i.e., with an NPT ensemble at 1.0 bar for 10 ns.In this process, the ethylene system was equilibrated at 90 K, in an amorphous solid state.Temperatures for anthracene are 140 and 260 K, and for pentancene, 50 and 300 K, set by the following experimental works. 49,50With all systems equilibrated in their solid states, crystalline and amorphous phases are both considered for anthracene and pentacene.The crystalline structures start from the corresponding X-ray structures. 71,72morphous structures were initialized using PACKMOL software (version 20.3.1). 73,74 RESULTS

Reorganization Energy from Ultraviolet
Photoelectron Spectroscopy.Ultraviolet photoelectron spectroscopy (UPS) has been employed to study the inner-sphere reorganization energy of anthracene and pentacene. 49,50With eqs 4 and 5, both low-frequency and high-frequency contributions of reorganization energy can be fitted from vibronic progression of UPS.Experimental spectra were fitted with a pseduo-Voigt function, which is a linear combination of Lorentzian and Gaussian functions, and fwhm values were reported.As shown in Table 1, λ s values range from 177 to 549 meV.The result from gas phase pentacene spectra was also included, and a much smaller λ s , 18 meV, from a much narrower width of vibronic bands was also included in Table 1.A gas-phase system does not have any outer-shell reorganization and the bandwidth serves as a lower bound for this approach.The width and the inferred λ s of solid-state UPS are well above this limit; therefore, inferring the broadening effect of the surrounding is not likely to be an instrumental artifact, and it should be reliable for inferring interactions with the environment.All λ s values are of the same order of magnitude as λ v , implying the importance of low-frequency reorganization energy, which is mainly from outer-sphere interactions.The temperature dependence of λ s is larger than that of λ v .Changing temperature could change the frequency of molecular vibrations, but its coupling to the charged state is not likely to change much.However, low-frequency modes contain intermolecular motions, which can be totally different at different temperatures; thus, their influence on the chargedstate energy can be quite significant.In this case we observed a reduction in higher temperature cases.

Ethylene at 90 K.
In Figure 2(a), the decay lifetime of ACF is fitted as 484 fs.As a comparison, the ACF for electronic coupling was reported in ref 75., and the corresponding decay lifetimes are 1083 fs (55K) and 210 fs (155 K).In Figure 2(b), the spectral density J(ω) of amorphous ethylene at 90 K is shown.There exist lowfrequency components (shown in the inset) together with several distinct peaks corresponding to their vibrational frequencies.λ out values calculated with eq 17 are listed in Table 2. λ out is obtained as 1088 meV.If eq 17 were integrated to 500 cm −1 , the low-frequency contribution is 968 meV (86%).This result implies that intermolecular motions are the dominant source of reorganization.

Anthracene at 140 and 260 K.
As seen in Figure 3, spectral densities are quite similar at different temperatures.Meanwhile, differences between amorphous and crystalline phases are more obvious, especially in the low-frequency region.The spectral density function is broader in the lowfrequency region in the amorphous phase, compared to those in the crystalline phase.
In Table 2, λ out calculated for crystalline and amorphous phases at 140 K are 212 and 358 meV, respectively, and lowfrequency contributions are 88 (42%) and 260 meV (73%).Since spectral density functions are similar between 140 and 260 K, reorganization energies at 260 K are also similar as at 140 K, which are 228 and 354 meV for crystalline and In units of meV.amorphous phases, respectively, and low-frequency contributions are 102 (45%) and 260 meV (73%).

Pentacene at 50 and 300 K.
In Figure 4, spectral densities J(ω) of crystalline and amorphous pentacene at 50K are included, again showing the importance of low-frequency modes.Like that for anthracene, the low frequency spectral density is broader in the amorphous phase than in the crystalline phase.
In Table 2, the λ out calculated for crystalline and amorphous phases at 50K are 170 and 184 meV, respectively, and the lowfrequency contribution is 84 (49%) and 108 meV (59%).At 300 K, values of λ out are 188 and 210 meV for crystalline and amorphous phases and the low-frequency contributions are 102 (54%) and 140 meV (67%).

DISCUSSION
4.1.General Values of λ out and λ s .Our results for λ s and λ out are summarized in Figure 5. Values of λ s , treated as an upper bound for λ out , are indeed higher than λ out calculated (the stacked height of the green bars) except for 300 K pentacene, where the two values are very close.We note that values of λ s derived from experimental spectra are from the crystalline phase.
The largest value of λ out , 1,088 meV, was for amorphous ethylene, the smallest molecule in the present work.Calculated λ out values for anthracene are all larger than that of pentacene.A λ out value that is higher than 1 eV maybe unusual.As the smallest model molecule with a π bond, such a value does not have any experimental implications either.Nevertheless, values of λ out , and λ s are all higher than 100 meV, which is about 4 times room temperature already, and they are certainly not negligible in practice.

High-Frequency Contribution of λ out in Nonpolar Systems.
A high-frequency contribution to λ out has rarely been discussed the literature, to the best of our knowledge.In traditional polar solvation works, vibrational peaks in infra absorption can be transformed into the dielectric function ε(ω) and peaks are seen in the spectral density function calculated from simple solvation models. 62,63However, in these cases, high-frequency contributions to the overall λ out are quite small, such that their quantum statistics do not really affect spectra of dynamics Stokes shift. 62For nonpolar systems, the same dielectric function yields very small λ out with narrow vibration peaks seen in the J(ω) modeled as well. 61Therefore, in the present work, λ out in the nonpolar systems simulated has a significant high-frequency component, as seen in Figure 5 and Table 2.This is a rather unexpected observation.Highfrequency components of λ out are similar, or even larger than low-frequency components in some cases, especially for the crystal phase (>50% for both temperatures in anthracene, and in pantenecene, 51% for 50K, and 46% for 300 K data.).
High-frequency components of λ out contribute to the third term on the right-handed side of Eq 8.The second term in Eq 8 can be estimated through density functional theory (DFT) for λ in with vibrational analysis, 76 and the values reported for both anthracene and pantacene are pretty small (0 or 1 meV for each mode, with 2 modes under 500 cm −1 ) in gas phase calculation, leading to a conclusion that real values of λ s can be much smaller than we have obtained from fitting UPS.However, MD and quantum chemistry are quite different in their model assumptions and sources of errors.A direct comparison and subtraction may not yield physically reliable values and further detailed studies may be necessary for a proper conclusion.Here we aim to obtain λ out and it is sufficient to conclude their significant magnitude.

Low Frequency Behavior of J(ω).
Environmental fluctuation to a system is often characterized using the Ohmic model for the spectral density. 77

J( ) exp( / )
s c s c 1 (19)   where ω c a characteristic "cut-off" frequency, and s indicates the power-law dependence at low frequencies.−80 Reporting fitted models from simulation offers practical model settings in these studies and better predictions and insights can be drawn.We fit our calculated spectral density function to Eq 19 and results are listed in Table 3.
As seen in Table 3, we found that all s values for the crystalline phase are larger than 1, indicating a super-Ohmic case, whereas s for amorphous systems are smaller than 1, Figure 4. Spectral density functions of pentacene in crystalline (Crys.)and amorphous (Amor.)phases.Blue and orange lines are for 50 and 300 K, respectively.Dark lines are results for the crystalline phase, and bright colors are for the amorphous phase.In the inset is a zoom-in for low frequency.Traces for the four data sets are shifted up for clarity.In the inset, their basal values are indicated by dashed lines of the same color, to the left.The vertical dashed line indicates the cutoff frequency at 200 cm −1 .
Figure 5. Experimental (Expt.) and calculated (Calc.)reorganization energies for both crystalline (Crys.)and amorphous (Amor.)states at the temperature indicated.λ out values are reported for calculated data, while for experimental data, inferred λ s values are included.Dark and light green colors indicate contributions to λ out from low-and highfrequency portions in spectral density, respectively.indicating a sub-Ohmic situation.As discussed earlier, in Figures 3 and 4, different behaviors in J(ω) in low-frequency limits are seen.The narrower low-frequency bands for J(ω) in crystallines are mainly in the left side of the band, where the slope in J(ω) goes from flat to steep before reaching the peak, showing a higher-order dependence of ω.
In ref 75, we reported the spectral density function for the fluctuation of electronic coupling, or the off-diagonal Hamiltonian matrix element.There we reported sub-Ohmic dynamics for amorphous phases.Our present findings of both sub-and super-Ohmic dynamics for fluctuations of systembath interactions in different phases can be useful for future studies on dynamics in ET systems.They also motivate us to reconsider sub-Ohmic dynamics we obtained for the offdiagonal Hamiltionian.It would be interesting to study fluctuation dynamics in electron transfer in a crystalline phase in the future.
In amorphous solid systems we can see that our MD simulated J(ω) does not approach zero at low frequency.In Figure 6, the lowest frequency data point for J(ω) is located at about 1.8 cm −1 and values of spectral densities from the correlation function of ΔE int , as well as those from rotational correlation, remain at a finite value in the amorphous cases.This is due to a long-time correlation from MD trajectories, especially derived from amorphous phases.This result indicates that the entire correlation of such movement is not fully captured and an ensemble of even longer trajectories is needed.However, for the purpose of the present work, in estimating λ out and further analyzing the physical origin, our present computational setting is sufficient.It would be interesting to further study such long-time correlation behavior in the future, since it allows us to further study properties related to inhomogeneous broadening, a set of important characteristics in amorphous solid states of material, but hard to study in detail.
4.4.Possible Sources of λ out .Our resulting λ out values are similar to those estimated for nonpolar solvents. 34,35The MD force field contains source physics of λ out .Since we are employing a classical force field, the effect of quadrupolar solvation is included. 35For the indirect effect of induced polarization changes with structure factor, 34 a polarizable force field would need to be employed, but it is not present in our current simulation, since ΔE int is simulated with only electrostatic interactions of the cation.Therefore, we believe that the observed λ out is from the quadrupolar or other higherorder moments.
Since the MD simulation was performed with the goal of determining λ out , arising from the electrostatic interaction between the cation and surrounding molecules, we did not include changes in structure when turning a neutral molecule into a cation, as they contribute to the inner-sphere reorganization energy.In so doing, we have ignored possible contributions to λ out from the movement of surrounding molecules in response to a change in molecular geometry of the cation.Without this contribution, albeit minor, λ out values we have obtained are not negligible already.
In addition, the previous study revealed the importance of rotational and translation motions in the low-frequency band of the spectral density of off-diagonal disorder. 75Therefore, dynamics, their corresponding ACF, and spectral densities were also investigated in the present work.To capture translational characteristics, the ACF is defined as where v⃗ is the velocity of central mass for each molecule.The rotational ACF, C rot (t), is defined through dynamics of the principle axis with the smallest moment of inertia, which is also aligned with the long molecular axis, as With P ⃗ being a normalized principle axis vector in the molecule.C rot (t) is an averaged cosine function of the angle between the two vectors.Figure 6 shows rotational and translation spectral densities, together with the spectral densities of ΔE int , and the fitted Ohmic function.The spectral density for ΔE int is a general combination of the two, with its higher frequency edge near 100 cm −1 following that of the translation spectral profile, and the lower frequency rotational.In the crystalline phase, all  three curves are very similar in the region of ∼1−30 cm −1 .In this region, slopes of these curves become steeper as frequency increases, indicating a clear super-Ohmic character.
It is somewhat surprising to see that rotational spectra in the amporphous phase have a significant contribution in a lower frequency region, as compared to translational spectra (1−10 cm −1 ), in all four cases studied.Very similar behavior for ΔE int is also seen in this region.This is consistent with data shown in Figure 5, where values of λ out are larger in the amorphous phase than the crystalline phase, mainly in their low-frequency contribution.Thus, it is likely that additional, low-frequency components in λ out are mainly constrained rotation.In an amorphous solid, without the ordered planar stacking in crystallines, molecules may have more space where hindered rotation is allowed.In this region, curves do not increase significantly before they reach the broad peak, indicating a sub-Ohmic characteristics.
4.5.Effects of λ out in Estimating ET Rates.ET rates with both the Marcus and Jortner's expressions (eqs 1 and 3) under different λ out (or λ s ), as a function of ΔE 0 are shown in Figure 7. λ out not only shifts the overall profile to a more negative ΔE 0 region, it also broadens the Gaussian profile for the Marcus rate.λ s has a similar effect for the Jortner's rate: the vibronic progression bands are moved toward a more negative ΔE 0 as they become broader.With zero λ out , it is not a problem to obtain a Marcus rate, but zero λ out implies a small λ s , and in the Jortner's expression, it would lead a sum of very narrow Gaussian functions.The steep dependence on ΔE 0 is not physically plausible, and it becomes hard to ensure predicted quality.A properly determined λ out (which is included in λ s ) is essential to predict ET rate, even for nonpolar systems.
4.6.Microscopic Dielectric Solvation.Our observation of significant outer-sphere organization energy in nonpolar systems is a deviation from bulk dielectric solvation models.A plausible way to understand this result is that, compared to the bulk, in a microscopic system, a larger dielectric solvation effect exists.This is consistent with many discussions for screening of charges in molecular dynamics, especially for biomolecules.For example, by studying the effect of polarizable solvation, dielectric constants that are larger than that from electronic polarization were reported. 81The dielectric constant formulation is related to the mean square of total dipole in a sampled volume, which contains the mean dipole, as well as the fluctuation of the dipole.

CONCLUSIONS
We studied outer-shell reorganization energies λ out in nonpolar molecular charge transporting systems and found that values can approach or even exceed 100 meV.Furthermore, the contribution at low-frequency is larger than 50%.This result indicates that translational and rotational modes of surrounding molecules are important in accommodating the charged molecule.Both crystalline and amorphous structures were also included in the simulation, and results show only a minor difference in the value of λ out .However, further analysis of their spectral density functions revealed sub-Ohmic dynamics for amorphous phase and super-Ohmic for crystalline phase.Nevertheless, this work demonstrates that λ out cannot be ignored, even in nonpolar systems.

Figure 1 .
Figure 1.(a) Depiction for different separations of λ: the traditional inner-and outer-shell separation in the right and left columns (λ in vs λ out ).The low-and high-frequency separations are depicted in upper and lower rows, which are modeled as λ s and λ v in fitting photoelectron spectra, as further schematically illustrated in (b).

Figure 2 .
Figure 2. (a) The autocorrelation function of the interaction energy ΔE int , C(t), for amorphous ethylene at 90 K; (b) the corresponding spectral density J(ω).In the inset is a zoom-in for low frequency.

Figure 3 .
Figure 3. Spectral density functions of anthracene in crystalline (Crys.)and amorphous (Amor.)phases.Blue and orange lines denote 140 and 260 K, respectively.Dark lines are results for the crystalline state, and bright colors are for the amorphous phase.In the inset is a zoom-in for low frequency.Traces for the four data sets are shifted up for clarity.In the inset, their basal values are indicated by dashed lines of the same color to the left.The vertical dashed line indicates the cutoff frequency at 200 cm −1 .

Figure 6 .
Figure 6.Spectral density curves of energetic disorder (black thick lines), their ohmic fitted curves (dashed lines), and those derived from translation (blue lines) and rotation (red lines) for (a) anthracene and (b) pentacene.In both (a) and (b), the upper and lower panels represent low and high temperatures, respectively, and the left and right columns depict crystalline and amorphous phases, respectively.

Figure 7 .
Figure 7. (a) ET rates with different λ out values (in eV) for the classical Marcus ET rate, eq 1, with a fixed λ in as 0.2 eV.(b) ET rates from the Jortner's expression, eq 3, with different λ s (in eV) as shown.λ v was set to 0.2 eV or with S being 1.169, and ℏω ν was set to 1380 cm −1 .Other parameters employed are V, 0.1 eV, and T, 50 K.

Table 1 .
High-and Low-Frequency Contribution of Reorganization Energies Inferred from UPS Spectra a

Table 2
a Calculated with trajectories of 18 ps, with each trajectory generating more than 300 sampled cations.b Contribution from frequencies under 500 cm −1 for ethylene or under 200 cm −1 for anthracene and pentacene.c In units of meV.