Breakdown of the Perturbative Approach to Molecular Packing

Jun 3, 2019 - Black and white dots represent the maxima of |VSF|2 and k, respectively. The peak positions of ∑m|Cmα|2|Cmβ|2 (white cross mark) and...
1 downloads 0 Views 4MB Size
Download PDF
Article Cite This: J. Phys. Chem. C 2019, 123, 15403−15411

pubs.acs.org/JPCC

Breakdown of the Perturbative Approach to Molecular Packing Dependence of Singlet Fission Rates in Pentacene Dimer Models: A Systematic Comparison with the Quantum Master Equation Approach Kenji Okada,† Takayoshi Tonami,† Takanori Nagami,† and Masayoshi Nakano*,†,‡,§,∥

Downloaded via UNIV PARIS-SUD on August 5, 2019 at 14:20:25 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

Department of Materials Engineering Science, Graduate School of Engineering Science, ‡Center for Spintronics Research Network (CSRN), Graduate School of Engineering Science, and §Quantum Information and Quantum Biology Division, Institute for Open and Transdisciplinary Research Initiatives, Osaka University, Toyonaka, Osaka 560-8531, Japan ∥ Institute for Molecular Science, 38 Nishigo-Naka, Myodaiji, Okazaki 444-8585, Japan S Supporting Information *

ABSTRACT: We have compared the molecular packing dependences of singlet fission (SF) rates calculated by perturbative and quantum master equation approaches using slip-stacked pentacene dimer models to clarify the applicability of the perturbative treatment and the importance of vibronic coupling effects on the SF dynamics. Several significant differences between these two approaches are found in the peak positions as well as in the relative peak amplitudes on the two-dimensional map of the SF rate as a function of the longitudinal and lateral displacements of the molecular backbones. Using the relative relaxation factor analysis, these disagreements are more striking for increasing the electronic couplings between the related diabatic exciton states as well as for reducing the energy difference between the lowest adiabatic Frenkel-like and double-triplet-like exciton states to match the primary vibronic coupling frequency. These results systematically illuminate the breakdown of the perturbative approach to the prediction of SF rates and indicate that exciton dynamics with explicit vibronic coupling is indispensable for obtaining a quantitative molecular packing dependence of SF rates. meV22−24) than the one-step mechanism by FE−TT (of the order of 1 meV22−24), although it cannot be a real pathway when the E(CT) is significantly larger than E(FE) and E(TT).22 It is predicted that when CT states lie higher in energy than the FE and TT states, but not excessively so, SF tends to occur via a superexchange mechanism virtually mediated by CT states,23,24 like in crystalline pentacene.25,26 Hence, the (Marcus-like) perturbation theory with the superexchange mechanism has often been applied to predicting SF rates27 and has been found to reproduce the trend of several experimental results,28 although it is known that there exist some limits of this kinetic model.29,30 Although nonperturbative approaches for exciton dynamics in multimers with multiple exciton states have been recently investigated,23,24,31−44 the comparison of the results between the SF dynamics simulation and the perturbative calculation has not been investigated in detail. Especially, there have been no systematic and comprehensive investigations on the theoretical rate of CT-mediated SF in typical test sets such as pentacene dimer models with various molecular packings, in

1. INTRODUCTION Singlet fission (SF) is a photophysical process in which a singlet exciton state (S1S0, S0S1), or a Frenkel exciton (FE) state, is converted into a correlated triplet (double-triplet) exciton state1(T1T1) (hereafter referred to as TT), which is followed by decoupling and spatial separation into two triplet states.1−4 Although SF was originally reported about 50 years ago,5,6 the interest in SF was revived in 2006, when Hanna and Nozik proposed multiple exciton generation through SF as a way to break the Shockley−Queisser limit7 on the efficiency of single-junction photovoltaics.8 Indeed, the explorations for novel SF materials as well as their mechanisms have been conducted intensely both experimentally and theoretically.1−4,9−20 Recently, various theoretical models for investigating the SF rate and TT yield have been developed in fields ranging from the energetics at the monomer level to the SF dynamics at the molecular aggregate level.2−4 One of the simplest approaches for this is to apply the perturbation theory (Fermi’s golden rule) (PT) or Marcus theory21 to evaluate the SF rate, which is composed of a one-step process (direct transition from FE to TT) or a two-step process via charge transfer (CT) states, that is, CA and AC (C: cation; A: anion). The two-step mechanism has larger effective couplings (of the order of 101−102 © 2019 American Chemical Society

Received: February 22, 2019 Revised: May 20, 2019 Published: June 3, 2019 15403

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C

Adiabatic exciton states {|α⟩} = {|FE′+⟩, |FE′−⟩, |CT′+⟩, |CT′−⟩, |TT′⟩} are described by the superposition of diabatic exciton states |FE+⟩ = √2/2(|S1S0⟩ + |S0S1⟩), |FE−⟩ = √2/2(| S1S0⟩ − |S0S1⟩), |CT+⟩ = √2/2(|CA⟩ + |AC⟩), |CT−⟩ = √2/ 2(|CA⟩ − |AC⟩) and |TT⟩, where the prime implies that the adiabatic state primarily includes the corresponding diabatic exciton configuration. 2.2. Relative Relaxation Factor Analysis. To clarify the exciton relaxation between the adiabatic exciton states, we employ the relative relaxation factor (RRF)41,43,46,47 in the case of Holstein coupling under the Markovian limit.

which completely the same parameters are employed in both perturbative and nonperturbative methods. That is why the scope of application of the traditional perturbative methods for clarifying the molecular packing dependences of SF rates has remained unclear and such perturbative methods have often been employed despite their being out of scope. Thus, as a suitable way to remedy this situation, we have planned a systematic comparison study, where the perturbative results must be evenly compared with the more accurate nonperturbative ones. In this study, therefore, we investigate the SF rates calculated by the perturbative method (Marcus-like approach) with the superexchange mechanism as well as by the time-convolutionless quantum master equation (QME) method41,43 to determine the reliability of the perturbative approach to the molecular configuration dependences of SF rates in slip-stacked pentacene dimer models, which are known to be typical SF systems.

ΔΓ(α → β) =

m

∑ γm(ω)

jijA (ω)ρ (t )A (ω)† − 1 {A (ω)† A (ω), ρ (t )}zyz m m m j m z ex ex 2 k { ω,m

Figure 1. Spectral density (J(|ωαβ|)) and relaxation parameters (γ(ωαβ), γ(ωβα), and γ(ωαβ) − γ(ωβα)) for T = 300 K and (Ω, λ) = (180, 50) meV. The net relaxation parameter γ(ωαβ) − γ(ωβα) has an extremum at |ωαβ| = Ω because it is proportional to the spectral density J(|ωαβ|).

(1)

Here, Am(ω) = ∑E(β)−E(α)=ω|α⟩⟨α|m⟩⟨m|β⟩⟨β|, where |α⟩ = ∑mCmα|m⟩ is an adiabatic state (an eigenstate of Hex with the eigenenergy E(α)) given by a linear combination of diabatic states {|m⟩} = {|S1S0⟩, |S0S1⟩, |CA⟩, |AC⟩, |TT⟩}. The second term on the right-hand side of eq 1 indicates the relaxation of the exciton density matrix ρex(t), which is characterized by the relaxation parameter γm(ω) given by the following equations (under the Markov approximation)41,43

adiabatic exciton states α and β include common diabatic exciton configurations m (|Cmα|2|Cmβ|2 ≠ 0), ΔΓ(α → β) has a large amplitude. In other words, the product of the vibronic coupling (spectral density) and the overlap of common diabatic exciton configurations between adiabatic exciton states determines the RRFs, which characterize SF dynamics. 2.3. Perturbation Method. The SF rate derived from the Marcus-like perturbation theory is determined by the square of the effective SF coupling, |VSF|2, at an arbitrary temperature under the assumption that the reorganization energies λ and the Gibbs free energy gap ΔG between S1S0 and TT are independent of the molecular configuration in the dimer. In this study, we employ diabatic exciton state energies obtained from the monomer calculation, so that λ and ΔG are regarded as substance-specific constants. This approach to the estimation of the SF rate has often been adopted in previous studies.27−29 Therefore, the relative SF rates for different molecular configurations can be discussed by comparing their | VSF|2. From the quasidegenerate perturbation theory, VSF is approximately expressed by the following equation2,48,49

γm(ω) = 2πJm (ω)(1 + nB(ω , T )) for ω > 0 λm kBT Ωm

for ω = 0

γm(ω) = 2πJm ( −ω)nB( −ω , T )

for ω < 0

γm(ω) = 4

(2)

where the spectral density Jm(ω) is of the Holstein vibrational mode of the mth diabatic exciton state and nB(ω, T) indicates the Bose−Einstein distribution at temperature T. The spectral density Jm(ω) is approximated by an Ohmic function with an identical Lorentz−Drude cutoff Ωm = Ω = 180 meV and reorganization energy λm = λ = 50 meV (typical values for the C−C stretching mode for acene molecules),23 defined in ω ≥ 0 with a maximum at ω = Ω Jm (ω) =

2 λmΩmω π ω 2 + (Ωm)2

(4)

where ωαβ = E(α) − E(β). The positive ΔΓ(α → β) value indicates the population damping from α to β, whereas the negative value indicates the population feeding from α to β. Although the relaxation parameter γm(ωαβ) includes the thermal distribution of the phonon, the net relaxation parameter γm(ωαβ) − γm(ωβα) is proportional to Jm(|ωαβ|). Thus, the amplitude of ΔΓ(α → β) takes a peak at |ωαβ| = Ω (see Figure 1). In addition, as seen from eq 4, when both the

2. METHODOLOGY 2.1. Quantum Master Equation Method. We consider the total Hamiltonian of the molecular dimer given by H = Hex + Hph + Hex−ph, where Hex is the exciton (electron) Hamiltonian; Hph and Hex−ph represent the phonon Hamiltonian and the exciton−phonon (vibronic) coupling Hamiltonian, respectively (see the Supporting Information). In this study, we consider only Holstein coupling, which is predicted to have significant effects on the SF dynamics.23,31−33 The second-order time-convolutionless master equation under the Markov approximation is expressed in atomic units (ℏ = m = e = 1) by the following equation41,43,45 d ρ (t ) = −i[Hex , ρex (t )] + dt ex

∑ |Cmα|2 |Cmβ|2 [γm(ωαβ) − γm(ωβα)]

(3) 15404

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C VSF =

1 2

ij yz 3 1 1 zz (VllVlh − VhhVhl)jjjj + zz 2 E (S S ) E (CT) E (TT) E (CT) − − 1 0 k {

3. RESULTS AND DISCUSSION 3.1. Molecular Packing Dependences of Singlet Fission by the Perturbation Theory and QME Dynamics. The square of the effective SF coupling |VSF|2 and the SF rate k on the ΔL−ΔT plane are shown in Figure 3a,b, respectively.

(5)

where Vij indicates the Fock matrix element between the molecular orbital (MO) i in molecule A and MO j in molecule B (i, j = h (the highest MO (HOMO) of the molecule), l (the lowest MO (LUMO) of the molecule)) (see the Supporting Information). 2.4. Pentacene Dimer Model. As shown in Figure 2, we consider slip-stacked pentacene dimer models in which the two

Figure 2. Top view (a) and side view (b) of the slip-stacked pentacene dimer model (0.00 Å ≤ ΔL ≤ 10.00 Å and 0.00 Å ≤ ΔT ≤ 2.75 Å). Coordinate axes are also shown.

Figure 3. 2D maps of squared effective SF couplings |VSF|2 [meV2] (a) and SF rates k [ps−1] (b) on the ΔL−ΔT plane. FE′− is found to be lower than TT′ in energy at the shaded area in (b). The white circle with a number denotes the presence of a peak. Peaks 1 and 2 in (b) are chosen to satisfy E(FE′−) > E(TT′), although they may not be the maxima.

molecules are stacked parallel to each other with a 3.50 Å separation (along the y axis) and the center of mass is shifted by ΔL and ΔT along the longitudinal (z) and transverse (x) axes, respectively, whose monomer structure is optimized by the spin-flip time-dependent density functional (SF-TD-DFT) method using the BHHLYP functional and 6-311G* basis set with Tamm−Dancoff and noncollinear approximations.50 The square of the effective coupling |VSF|2 and the SF rate k are investigated within the range of 0.00 Å ≤ ΔL ≤ 10.00 Å and 0.00 Å ≤ ΔT ≤ 2.75 Å, whereas electronic couplings and other values are calculated within ΔL ≤ 12.00 Å. The effective coupling and the matrix elements of the exciton Hamiltonian (energies of diabatic exciton states (diagonal elements) and electronic couplings (off-diagonal elements)) are obtained by the DFT and TD-DFT levels of approximation with the longrange corrected (LC)-BLYP functional (the range separating the parameter μ is set to 0.20 bohr−1 as determined by the IP tuning method51) and the 6-311G* basis set. For the diagonal elements, we assume that E(S1S0) (= E(S0S1)) is equal to the optically allowed first excitation energy E(S1) = 2260 meV and the TT state excitation energy can be approximated by E(TT) ≈ 2E(T1 ) = 1720 meV. These values are found to semiquantitatively reproduce the values employed in the previous SF dynamics study.23 The CT state excitation energy E(CT) (≡ E(CA) = E(AC)) is approximately estimated by eq S12 in the Supporting Information.34 These quantum chemical calculations are performed by the Gaussian 09 package.52 The TT population is obtained by applying the six-order Runge−Kutta method to solve the QME (eq 1) with the initial state S1S0 at T = 300 K. The SF rate k is estimated by fitting the time-dependent TT population PTT(t) with a threeparameter function PTT(t) = PTT(∞) − A exp(−kt), where PTT(∞) corresponds to the TT yield at infinite time. Note here that A = PTT(∞) − PTT(0) is satisfied within the numerical error since the initial population of S1S0 is 1.0 in the present case.

The molecular packing dependence of the singlet fission rate is found to be in semiquantitative agreement with the previous result.34 It is found that in the shaded area in Figure 3b, the adiabatic state with the largest population at infinitely long time is not a TT-like but an FE-like state primarily including S1S0 and S0S1 contributions. The trend of the two-dimensional (2D) map of |VSF|2 is found to support the SF rate k obtained by the SF dynamics simulation using the QME method. However, non-negligible differences are found in (i) the peak positions and (ii) the relative peak amplitudes. For example, as seen from Table 1, peak 3(Q) appears 0.5 Å distant in the ΔT direction from the corresponding 3(P) position. Moreover, peak 1(Q) gives k = 18 ps−1, which is 3.1 times as large as that at 1(P), in spite of a much smaller effective coupling at 1(Q) (only 14% of the value at 1(P)). Likewise, k = 22 ps−1 is found at 2(Q) (3.0 times as fast as that at 2(P)) despite the smaller |VSF|2 (see Table 1). These disagreements of the squared effective coupling and the SF rate are observed particularly when |VSF| is on the order of 10 meV or larger, which is reported to make the Marcus-like perturbation theory invalid.23,29 In addition to these strongly coupled regions, unfortunately, it is found to be difficult to obtain quantitative agreement between the perturbation theory and QME dynamics even in weakly coupled regions (less than 10−20 meV of |VSF|) (see the Supporting Information for details). Also, in the shaded region in Figure 3b with ΔL ∼ 2 Å, ultrafast exciton dynamics is predicted despite the small |VSF| values (less than the order of 10 meV). The large variations of the ratios k/|VSF|2 at different peaks listed in Table 1 indicate that the SF rate k is not usually proportional to |VSF|2 but that the ratios strongly depend on the peak position. This large variation feature in k/|VSF|2 is shown to be more distinct for smaller ΔL values, whereas 15405

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C

Table 1. SF Rates, Effective SF Couplings, and TT Yield at the Peak Position (ΔL, ΔT) Evaluated by the Perturbation Theory (P) and QME Dynamics (Q) peaka 1 2 3 4 5 6

(P) (Q)b (P) (Q)b (P) (Q) (P) (Q) (P), (Q)c (P), (Q)c

ΔL [Å]

ΔT [Å]

k [ps−1]

|VSF|2 [×103 meV2]

k/|VSF|2 [×10−4 ps−1 meV−2]

PTT(∞) [−]

2.00 2.25 3.25 3.00 4.50 4.50 5.75 5.50 6.75 9.00

1.50 0.75 0.75 0.50 1.50 1.00 0.75 0.50 2.00 2.00

5.85 18 7.44 22 2.95 3.44 0.378 0.510 0.845 0.0871

4.95 0.680 3.44 2.89 3.67 2.82 0.458 0.219 1.03 0.107

1.18 26 2.16 7.6 0.804 1.22 0.820 2.33 0.820 0.814

0.933 0.699 0.942 0.685 0.936 0.954 0.990 0.995 0.967 0.994

a

These numbers correspond to those in Figure 3. b(Q) is chosen to satisfy E(FE′−) > E(TT′). cBoth peaks (P) and (Q) coincide with each other.

k/|VSF|2 values turn out to converge for larger ΔL values, where the effective coupling is shown to be sufficiently small. As a result, the positions of peaks 5(P) and 6(P) on the |VSF|2 map are found to coincide with those of the corresponding (Q) on the k map. These findings indicate that the perturbative approach is relatively reliable only in the region where all couplings are smaller than the order of 10 meV,23,29 but such a region is not important in the SF application due to its very small SF rate. 3.2. Relative Relaxation Factor Analysis of QME Results. We investigate the origin of the disagreement between these two methods (P) and (Q) by the RRF analysis. As shown in eq 4, the product of the diabatic exciton overlap between two adiabatic exciton states ∑m|Cmα|2|Cmβ|2 and the relaxation parameter γ (related to the vibronic coupling) defines the RRFs. Here, the lower FE-like adiabatic state (FE′+ or FE′−) and the TT′ state (see also Section 3.3) are indicated by α and β, respectively. The SF dynamics is known to be approximately described by the population relaxation between these two adiabatic states.39,41 Indeed, the map of ΔΓ(α → β) (Figure 4a) is shown to be similar to that of k (Figure 3b) in the peak positions and the relative peak amplitudes. ΔΓ(α → β) shows negative values for ΔL < 1.0 Å and ΔT < 1.0 Å (shown in purple) because α (FE′−) is found to be lower than β (TT′) in that region. In this case, the relaxation from FE′+ to FE′− is found to be dominant since both the states include common FE (S1S0, S0S1) configurations. Figure 4b,c show the overlap of diabatic exciton configurations ∑m|Cmα|2|Cmβ|2 and the net relaxation parameter γ(ωαβ) − γ(ωβα). Figure 4b indicates that contributions of the diabatic exciton overlap determine the approximate landform of the 2D map of SF rates, while peak positions of ∑m|Cmα|2|Cmβ|2 are different from those of the perturbation theory for ΔL ≤ 6.0 Å. Concerning peak 3, for example, the maximum is shown to be located between 3(P) and 3(Q). In addition, the exciton overlap between α and β is found to be primarily caused by TT and FE configurations when the state α is FE′−. These features originate in the perturbative treatment of electronic coupling in the region where the approach of FE′− and TT′ energies enhances both TT configuration in FE′− and FE configurations in TT′ (see Figures S8 and S9 in the Supporting Information), although these deficiencies in the perturbative treatment may be improved by using the higher-order perturbation theory. The increase of γ(ωαβ) − γ(ωβα) value with decreasing ΔT in the (ΔT ≤ 2.0 Å, ΔL ≤ 6.0 Å) region except for (ΔL < 1.0 Å, ΔT < 1.2 Å) and (2.0 Å < ΔL < 3.2 Å, ΔT < 1.0 Å) indicates the

Figure 4. Relative relaxation factor ΔΓ(α → β) (a), the overlap between α and β (b), and the net relaxation parameter γ(ωαβ) − γ(ωβα) (c). Black and white dots represent the maxima of |VSF|2 and k, respectively. The peak positions of ∑m|Cmα|2|Cmβ|2 (white cross mark) and |VSF|2 do not coincide with each other for ΔL ≤ 6.0 Å; for example, the position giving the maximum of ∑m|Cmα|2|Cmβ|2 for peak 3 is located between 3(P) and 3(Q).

increase of the vibronic coupling contribution to the SF rate with decreasing ΔT, which is found to cause the relative increase of the SF rate by the QME method as compared with that by the perturbation theory. As shown in Table 2, ∑m|Cmα|2|Cmβ|2 becomes larger at (Q) than at (P) for peaks 1, 2, and 4, but 3(Q) has a slightly smaller ∑m|Cmα|2|Cmβ|2 value. On the other hand, γ(ωαβ) − γ(ωβα) value at peaks 3(Q) and 4(Q) are found to be larger than their counterparts at 3(P) and 4(P), respectively. These results indicate that (i) the significantly large exciton overlap contribution causes the shift of peaks 1 and 2, (ii) the vibronic coupling causes the shift of peak 3, and (iii) both the exciton overlap and the vibronic coupling cause the shift of peak 4. It is also interesting that FE′+ is found to be the state α at 4(P), 15406

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C

Table 2. Relative Relaxation Factors and Their Components at the Peak Position (ΔL, ΔT) Evaluated by the Perturbation Theory (P) and QME Dynamics (Q) peaka 1 2 3 4 5 6

(P) (Q)b (P) (Q)b (P) (Q) (P) (Q) (P), (Q)c (P), (Q)c

ΔL [Å]

ΔT [Å]

ΔΓ(α → β) [meV]

∑m|Cmα|2|Cmβ|2 [−]

γ(ωαβ) − γ(ωβα) [meV]

2.00 2.25 3.25 3.00 4.50 4.50 5.75 5.50 6.75 9.00

1.50 0.75 0.75 0.50 1.50 1.00 0.75 0.50 2.00 2.00

4.04 12 5.36 23 2.16 2.36 0.105 0.339 0.645 0.0597

0.0546 0.228 0.0601 0.271 0.0342 0.0334 0.00149 0.00408 0.0107 0.000977

74.0 51 89.2 83 63.3 70.6 70.4 83.1 60.4 61.1

a

These numbers correspond to those in Figure 3. b(Q) is chosen to satisfy E(FE′−) > E(TT′). cBoth peaks (P) and (Q) coincide with each other.

Figure 5. Energy diagrams of adiabatic exciton states and electronic couplings at peak 1(P) (ΔL, ΔT) = (2.00, 1.50) Å (a) and 1(Q) (ΔL, ΔT) = (2.25, 0.75) Å (b).

low-frequency modes of intermolecular Peierls coupling, which has often been reported to have significant effects on the SF rate36,54−56 although it has been ignored in this study, will be one of the issues to be considered in the future. It should be reemphasized that the remarkable increase in SF rates is observed around 1(Q) and 2(Q), which is not expected from the perturbation theory. This marked phenomenon is found to originate from the small α−β energy gap, which causes enormous mixing of diabatic configurations in α and β, although the energy gap is sometimes found to be too small to gain the benefit of the vibronic coupling contribution. At the same time, however, a remarkable decrease in the TT yield is expected in such a case (see Table 1). This result is found to coincide with the tendency that (nearly) perfect energy matching between diabatic FE and TT states enhances the SF rates but significantly lowers TT yields, which was found in our previous study.41 3.3. Analysis of Molecular Packing Dependences of Adiabatic Exciton States and Diabatic Exciton Configurations Based on Electronic Couplings. The diabatic exciton configurations in the adiabatic exciton states and the energies of the adiabatic exciton states discussed above are described by electronic couplings. We consider an alternative diabatic basis set {|FE+⟩ = √2/2(|S1S0⟩ + |S0S1⟩), |FE−⟩ = √2/2(|S1S0⟩ − |S0S1⟩), |CT+⟩ = √2/2(|CA⟩ + |AC⟩), |CT−⟩ = √2/2(|CA⟩ − |AC⟩), |TT⟩}, which can diagonalize each diabatic exciton block FE (S1S0 and S0S1), CT (CA and AC), and TT. The Hex in this alternative basis set representation is expressed as follows

which causes a small exciton overlap between FE′+ and TT′ because FE′+ does not include the diabatic TT configuration. We next investigate the origin of the difference in the relative amplitude of each peak by focusing on peaks 2 and 3. As seen from Table 1, almost the same |VSF|2 values are obtained at 2(P) and 3(P), and also 2(Q) and 3(Q), whereas the QME approach predicts significantly faster dynamics at 2(Q) and 2(P) than at 3(Q) and 3(P), respectively. As seen from Table 2, the RRF at 2(Q) is 9.5 times as large as that at 3(Q), where the RRF at 2(Q) is composed of 8.1 times larger ∑m|Cmα|2|Cmβ|2 and 1.2 times larger γ(ωαβ) − γ(ωβα) values than those at 3(Q), respectively. It is similarly found that the RRF at 2(P) is 2.5 times larger than that at 3(P) because of the contributions of the 1.8 times larger exciton overlap and 1.4 times larger vibronic coupling at 2(P) than at 3(P). These features are understood by the fact that the electronic coupling that lowers E(FE′−) and increases TT configuration in FE′− tends to be larger for smaller ΔL values. Therefore, the synergistic effect of exciton overlap (configuration mixing) and the vibronic coupling is found to have a great impact on the relative rate amplitude of the two maxima of peaks 2 and 3. Likewise, the differences in the relative amplitude between 1 and 3 are found to be described using these two components of the RRFs. Since the state-dependent vibronic coupling has been reported to affect the SF dynamics53 (although we ignore the effect in this study for simplicity), the present finding on the relationship between the vibronic coupling and SF dynamics is predicted to be more important when different spectral densities are given for each diabatic state. 41 Furthermore, strong vibronic coupling, particularly in the 15407

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C ij E(S1) + Vex jj jj jj 0 jj jj j Hex = jjjjVll − Vhh jj jj jj 0 jj jj jj k0

0 E(S1) − Vex 0 Vll + Vhh 0

yz zz zz zz 0 Vll + Vhh 0 zz zz z E(CT) 0 3 /2(Vlh + Vhl)zzzz zz zz 0 E(CT) 3 /2(Vlh − Vhl)zzz zz zz z 3 /2(Vlh + Vhl) 3 /2(Vlh − Vhl) E(TT) {

Vll − Vhh

where Vex is defined by eq S9 and causes Davydov splitting. In the present case, √3/2(Vlh − Vhl) = 0 meV and √3/2(Vlh + Vhl) = √3Vlh. As shown in Figure 5, Vll − Vhh can affect CT′+ and FE′+ energies and their diabatic exciton configurations. Likewise, Vll + Vhh, and √3/2(Vlh − Vhl) are found to affect configurations of FE′+, CT′−, and TT′. We found that Vll + Vhh = −284.6 + (−70.2) = −354.8 meV at 1(P), whereas Vll + Vhh = −356.6 + (−347.6) = −704.2 meV at 1(Q). Especially, the significant increase of |Vll + Vhh| value at 1(Q) is found to stabilize FE′−, but destabilizes CT′−. Likewise, the mixing and energy of CT′− and TT′ are affected by √3/2(Vlh − Vhl), which is found to have a relatively small amplitude: −207.8 meV at 1(P) and −37.3 meV at 1(Q). Besides, the energy gap between diabatic CT− and TT states is found to be larger than the FE−−TT gap. Thus, the influence of √3/2(Vlh − Vhl) on TT′ turns out to be modest as compared with the case of Vll + Vhh on FE′− (see Figure 5). However, it is notable that a nonzero √3/2(Vlh − Vhl) value can cause the mixing of FE to TT′ and, simultaneously, of TT to FE′− when the FE and TT energy levels are close to each other. Such a CT-mediated FE/ TT configuration mixing is not described by the conventional perturbation theory (eq 5). We have shown 2D maps of off diagonal elements of eq 6 in Figure 6a−c and those of FE, CT, and TT characteristics of adiabatic states (FE′− and TT′) in Figure 7. The distribution of the CT characteristic in FE′− on the ΔL−ΔT plane is

0

0

(6)

similar to the distribution of |Vll + Vhh| (Figure 6b). Likewise, the distribution of the CT characteristic in TT′ and that of √3/2|Vlh − Vhl| (Figure 6c) are found to be alike. It is found that the mixing of TT to FE′− and that of FE to TT′ tend to increase in the region where a larger |VSF| is found (Figures 3a and S5) because they are all affected by Vll + Vhh and √3/2(Vlh − Vhl) = √3Vlh (see eq 5 and Figure 5), although their mixing amplitudes are not completely the same. As mentioned above, inclusion of the TT configuration in FE′− and the FE configurations in TT′ cannot be described by the conventional perturbation theory (eq 5), whereas it is found that FE′− and TT′ obtained by diagonalization of the Hex matrix involve TT and FE configurations, respectively. As a consequence, among the exciton overlap contributions of each diabatic state to ΔΓ(FE′− → TT′), the contribution of the TT configuration is found to be the largest (Figure S9). For example, at (ΔL, ΔT) = (4.50, 1.25) Å, the TT contribution is found to be ∼5 and ∼15 times as large as the FE and CT contributions, respectively. Especially, mixing of the TT configuration to FE′− and the FE configuration to TT′ can be enhanced when FE′− and TT′ are near-degenerate. Because this configuration mixing cannot be described by the conventional perturbation theory (eq 5), it is predicted that the perturbative approach cannot reproduce SF rates obtained from a non-perturbative approach like the QME approach even semiquantitatively, particularly in strongly electronic-coupled systems, which are expected to be efficient SF systems. It is worthy to note that the electronic couplings, which affect energies and configurations of adiabatic states, can be understood by the phase and overlap of molecular orbitals. When recurring units in the HOMO of molecule A (distributions on the facing vertices of the upper and lower carbon zigzag edges, see Figure S1) are located just above those in molecule B, the constructive hA−hB overlap results in a large amplitude of Vhh (see Figure S4). In contrast, Vhh = 0 is found at ΔT ∼ 1.75 Å because the HOMO of pentacene has a nodal plane orthogonal to the transverse direction of the molecule (the x axis in Figure 2). Similarly, a large spatial overlap of lA−lB causes the large amplitude of Vll around every 2.25 Å to the longitudinal direction. As shown in Figure S1, because there are no nodal planes orthogonal to the transverse axis in the LUMO of pentacene, Vll turns out to increase its amplitude with decreasing ΔT. Due to the symmetry of the present dimer configuration, Vhl = −Vlh is found at any molecular configurations, simultaneously, and Vhl and Vlh are found to be zero at ΔT = 0 Å. When opposite zigzag edges of upper and lower pentacene molecules are overlapped, the amplitudes of Vhl and Vlh are found to increase because of the constructive-overlap contribution as well as of the smaller destructive contribution.

Figure 6. 2D maps of off-diagonal elements of eq 6: Vll − Vhh (a), Vll + Vhh (b), and √3/2(Vlh − Vhl) (c) on the ΔL−ΔT plane. The contour interval is set to 100 meV, and the contour of 0 meV is shown by bold lines in each figure.

4. CONCLUSIONS In this study, we have systematically illuminated the breakdown of the Marcus-like PT approach to predicting the 15408

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C

Figure 7. 2D maps of the FE, CT, and TT characteristics of adiabatic states (FE′− and TT′) on the ΔL−ΔT plane.

■

molecular packing dependence of SF rates in pentacene dimer models by comparing with the QME results. Completely the same parameters have been employed in both methods. Although the obtained 2D map of effective coupling |VSF|2(PT) is found to support the trend of the result of the SF rate k(QME) approximately, the maxima of |VSF|2 and k are found at different positions with different amplitudes. From the RRF analysis in the QME approach, we have discussed the origin of disagreements found by partitioning the relaxation rate into two factors: diabatic exciton overlap (mixing) between two adiabatic exciton states and the vibronic coupling effect. It is notable that when the electronic coupling between the diabatic states causes such crucial variations in the mixing of diabatic states, and a non-negligible difference exists between the peak frequency of the spectral density in vibronic coupling and the adiabatic exciton state energy gap, the conventional perturbative treatment has the possibility of failing to predict the SF rate quantitatively or even qualitatively. These problems tend to be inevitable for designing efficient SF materials with a high SF rate. In summary, the present results demonstrate that although the effective SF coupling obtained by the conventional perturbation theory can work as a simple indicator to just recognize the presence of the maxima of the SF rate in the search area of optimal molecular packing, a nonperturbative treatment like the QME approach with explicit vibronic coupling is indispensable for revealing the SF mechanism, predicting the quantitative SF rate and TT yield, and for designing novel SF materials.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Kenji Okada: 0000-0002-2185-0309 Takayoshi Tonami: 0000-0001-8518-9284 Takanori Nagami: 0000-0001-9252-0596 Masayoshi Nakano: 0000-0002-3544-1290 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by JSPS KAKENHI Grant Number JP18H01943 in Scientific Research (B), Grant Number JP17H05157 in Scientific Research on Innovative Areas “πSystem Figuration”, and Grant Number JP26107004 in Scientific Research on Innovative Areas “Photosynergetics”, and JSPS Research Fellowship for Young Scientists (JP18J20887 (T.N.)). Theoretical calculations were partly performed using Research Center for Computational Science, Okazaki, Japan.

■

REFERENCES

(1) Smith, M. B.; Michl, J. Singlet Fission. Chem. Rev. 2010, 110, 6891−6936. (2) Smith, M. B.; Michl, J. Recent Advances in Singlet Fission. Annu. Rev. Phys. Chem. 2013, 64, 361−386. (3) Casanova, D. Theoretical Modeling of Singlet Fission. Chem. Rev. 2018, 118, 7164−7207. (4) Japahuge, A.; Zeng, T. Theoretical Studies of Singlet Fission: Searching for Materials and Exploring Mechanisms. ChemPlusChem 2018, 83, 146. (5) Geacintov, N.; Pope, M.; Vogel, F. Effect on Magnetic Field on the Fluorescence of Tetracene Crystals: Exciton Fission. Phys. Rev. Lett. 1969, 22, 593−596. (6) Merrifield, R. E.; Avakian, P.; Groff, R. P. Fission of Singlet Excitons into Pairs of Triplet Excitons in Tetracene Crystals. Chem. Phys. Lett. 1969, 3, 386−388.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.9b01713. Total Hamiltonian for the pentacene dimer model; estimation of the SF rate from the perturbation theory; relaxation process in SF dynamics; contributions of respective diabatic exciton states to RRFs (PDF) 15409

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C (7) Shockley, W.; Queisser, H. Detailed Balance Limit of Efficiency of P-N Junction Solar Cells. J. Appl. Phys. 1961, 32, 510−519. (8) Hanna, M. C.; Nozik, A. J. Solar Conversion Efficiency of Photovoltaic and Photoelectrolysis Cells with Carrier Multiplication Absorbers. J. Appl. Phys. 2006, 100, No. 074510. (9) Paci, I.; Johnson, J. C.; Chen, X.; Rana, G.; Popović, D.; David, D. E.; Nozik, A. J.; Ratner, M. A.; Michl, J. Singlet Fission for DyeSensitized Solar Cells: Can a Suitable Sensitizer Be Found? J. Am. Chem. Soc. 2006, 128, 16546−16553. (10) Minami, T.; Nakano, M. Diradical Character View of Singlet Fission. J. Phys. Chem. Lett. 2012, 3, 145−150. (11) Minami, T.; Ito, S.; Nakano, M. Fundamental of DiradicalCharacter-Based Molecular Design for Singlet Fission. J. Phys. Chem. Lett. 2013, 4, 2133−2137. (12) Ito, S.; Minami, T.; Nakano, M. Diradical Character Based Design for Singlet Fission of Condensed-Ring Systems with 4nπ Electrons. J. Phys. Chem. C 2012, 116, 19729−19736. (13) Minami, T.; Ito, S.; Nakano, M. Theoretical Study of Singlet Fission in Oligorylenes. J. Phys. Chem. Lett. 2012, 3, 2719−2723. (14) Zeng, T.; Ananth, N.; Hoffmann, R. Seeking Small Molecules for Singlet Fission: A Heteroatom Substitution Strategy. J. Am. Chem. Soc. 2014, 136, 12638−12647. (15) Eaton, S. W.; Miller, S. A.; Margulies, E. A.; Shoer, L. E.; Schaller, R. D.; Wasielewski, M. R. Singlet Exciton Fission in Thin Films of tert-Butyl-Substituted Terrylenes. J. Phys. Chem. A 2015, 119, 4151−4161. (16) Ito, S.; Nagami, T.; Nakano, M. Diradical Character-Based Design for Singlet Fission of Bisanthene Derivatives: Aromatic-Ring Attachment and π-Plane Twisting. J. Phys. Chem. Lett. 2016, 7, 3925− 3930. (17) López-Carballeira, D.; Casanova, D.; Ruipérez, F. Theoretical Design of Conjugated Diradicaloids as Singlet Fission Sensitizers: Quinones and Methylene Derivatives. Phys. Chem. Chem. Phys. 2017, 19, 30227−30238. (18) Shen, L.; Wang, X.; Liu, H.; Li, X. Tuning the Singlet Fission Relevant Energetic Levels of Quinoidal Bithiophene Compounds by Means of Backbone Modifications and Functional Group Introduction. Phys. Chem. Chem. Phys. 2018, 20, 5795−5802. (19) Messelberger, J.; Grünwald, A.; Pinter, P.; Hansmann, M. M.; Munz, D. Carbene Derived Diradicaloids − Building Blocks for Singlet Fission? Chem. Sci. 2018, 9, 6107−6117. (20) López-Carballeira, D.; Casanova, D.; Ruipérez, F. Potential Use of Squarates and Croconates as Singlet Fission Sensitizers. Chem. Phys. Chem. 2018, 19, 2224−2233. (21) Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I. J. Chem. Phys. 1956, 24, 966−978. (22) Renaud, N.; Sherratt, P. A.; Ratner, M. A. Mapping the Relation between Stacking Geometries and Singlet Fission Yield in a Class of Organic Crystals. J. Phys. Chem. Lett. 2013, 4, 1065−1069. (23) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. Microscopic Theory of Singlet Exciton Fission. II. Application to Pentacene Dimers and the Role of Superexchange. J. Chem. Phys. 2013, 138, No. 114103. (24) Mirjani, F.; Renaud, N.; Gorczak, N.; Grozema, F. C. Theoretical Investigation of Singlet Fission in Molecular Dimers: The Role of Charge Transfer States and Quantum Interference. J. Phys. Chem. C 2014, 118, 14192−14199. (25) Sebastian, L.; Weiser, G.; Bässler, H. Charge Transfer Transitions in Solid Tetracene and Pentacene Studied by Electroabsorption. Chem. Phys. 1981, 61, 125−135. (26) Yamagata, H.; Norton, J.; Hontz, E.; Olivier, Y.; Beljonne, D.; Brédas, J. L.; Silbey, R. J.; Spano, F. C. The Nature of Singlet Excitons in Oligoacene Molecular Crystals. J. Chem. Phys. 2011, 134, No. 204703. (27) Bhattacharyya, K.; Datta, A. Polymorphism Controlled Singlet Fission in TIPS-Anthracene: Role of Stacking Orientation. J. Phys. Chem. C 2017, 121, 1412−1420.

(28) Bae, Y. J.; Kang, G.; Malliakas, C. D.; Nelson, J. N.; Zhou, J.; Young, R. M.; Wu, Y.-L.; Van Duyne, R. P.; Schatz, G. C.; Wasielewski, M. R. Singlet Fission in 9,10-Bis(phenylethynyl)anthracene Thin Films. J. Am. Chem. Soc. 2018, 140, 15140−15144. (29) Yost, S. R.; Lee, J.; Wilson, M. W. B.; Wu, T.; McMahon, D. P.; Parkhurst, R. R.; Thompson, N. J.; Congreve, D. N.; Rao, A.; Johnson, K.; et al. A Transferable Model for Singlet-Fission Kinetics. Nat. Chem. 2014, 6, 492−497. (30) Le, A. K.; Bender, J. A.; Arias, D. H.; Cotton, D. E.; Johnson, J. C.; Roberts, S. T. Singlet Fission Involves an Interplay between Energetic Driving Force and Electronic Coupling in Perylenediimide Films. J. Am. Chem. Soc. 2018, 140, 814−826. (31) Greyson, E. C.; Vura-Weis, J.; Michl, J.; Ratner, M. A. Maximizing Singlet Fission in Organic Dimers: Theoretical Investigation of Triplet Yield in the Regime of Localized Excitation and Fast Coherent Electron Transfer. J. Phys. Chem. B 2010, 114, 14168−14177. (32) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. Microscopic Theory of Singlet Exciton Fission. I. General Formulation. J. Chem. Phys. 2013, 138, No. 114102. (33) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. Microscopic Theory of Singlet Exciton Fission. III. Crystalline Pentacene. J. Chem. Phys. 2014, 141, No. 074705. (34) Wang, L.; Olivier, Y.; Prezhdo, O. V.; Beljonne, D. Maximizing Singlet Fission by Intermolecular Packing. J. Phys. Chem. Lett. 2014, 5, 3345−3353. (35) Tao, G. Bath Effect in Singlet Fission Dynamics. J. Phys. Chem. C 2014, 118, 27258−27264. (36) Renaud, N.; Grozema, F. C. Intermolecular Vibrational Modes Speed up Singlet Fission in Perylenediimide Crystals. J. Phys. Chem. Lett. 2015, 6, 360−365. (37) Tamura, H.; Huix-Rotllant, M.; Burghardt, I.; Olivier, Y.; Beljonne, D. First-Principles Quantum Dynamics of Singlet Fission: Coherent Versus Thermally Activated Mechanisms Governed by Molecular π Stacking. Phys. Rev. Lett. 2015, 115, No. 107401. (38) Zeng, T.; Goel, P. Design of Small Intramolecular Singlet Fission Chromophores: An Azaborine Candidate and General Small Size Effects. J. Phys. Chem. Lett. 2016, 7, 1351−1358. (39) Zeng, T. Through-Linker Intramolecular Singlet Fission: General Mechanism and Designing Small Chromophores. J. Phys. Chem. Lett. 2016, 7, 4405−4412. (40) Zheng, J.; Xie, Y.; Jiang, S.; Lan, Z. Ultrafast Nonadiabatic Dynamics of Singlet Fission: Quantum Dynamics with the Multilayer Multiconfigurational Time-Dependent Hartree (Ml-Mctdh) Method. J. Phys. Chem. C 2016, 120, 1375−1389. (41) Nakano, M.; Ito, S.; Nagami, T.; Kitagawa, Y.; Kubo, T. Quantum Master Equation Approach to Singlet Fission Dynamics of Realistic/Artificial Pentacene Dimer Models: Relative Relaxation Factor Analysis. J. Phys. Chem. C 2016, 120, 22803−22815. (42) Fujihashi, Y.; Chen, L.; Ishizaki, A.; Wang, J.; Zhao, Y. Effect of High-Frequency Modes on Singlet Fission Dynamics. J. Chem. Phys. 2017, 146, No. 044101. (43) Nakano, M.; Nagami, T.; Tonami, T.; Okada, K.; Ito, S.; Kishi, R.; Kitagawa, Y.; Kubo, T. Quantum Master Equation Approach to Singlet Fission Dynamics in Pentacene Linear Aggregate Models: Size Dependences of Excitonic Coupling Effects. J. Comput. Chem. 2019, 40, 89−104. (44) Nakano, M.; Okada, K.; Nagami, T.; Tonami, T.; Kishi, R.; Kitagawa, Y. Monte Carlo Wavefunction Approach to Singlet Fission Dynamics of Molecular Aggregates. Molecules 2019, 24, 541. (45) Breuer, H.-P.; Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press: Oxford, U.K., 2002. (46) Takahata, M.; Nakano, M.; Fujita, H.; Yamaguchi, K. Mechanism of Exciton Migration of Dendritic Molecular Aggregate: A Master Equation Approach Including Weak Exciton-Phonon Coupling. Chem. Phys. Lett. 2002, 363, 422−428. (47) Nakano, M.; Takahata, M.; Yamada, S.; Yamaguchi, K.; Kishi, R.; Nitta, T. Exciton Migration Dynamics in a Dendritic Molecule: Quantum Master Equation Approach Using Ab Initio Molecular 15410

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411

Article

The Journal of Physical Chemistry C Orbital Configuration Interaction Method. J. Chem. Phys. 2004, 120, 2359−2367. (48) Shavitt, I.; Redmon, L. T. Quasidegenerate perturbation theories. A Canonical van Vleck Formalism and Its Relationship to Other Approaches. J. Chem. Phys. 1980, 73, 5711−5717. (49) Ito, S.; Nagami, T.; Nakano, M. Design Principles of Electronic Couplings for Intramolecular Singlet Fission in Covalently-Linked Systems. J. Phys. Chem. A 2016, 120, 6236−6241. (50) Minami, T.; Ito, S.; Nakano, M. Signature of Singlet Open-Shell Character on the Optically Allowed Singlet Excitation Energy and Singlet−Triplet Energy Gap. J. Phys. Chem. A 2013, 117, 2000−2006. (51) Baer, R.; Livshits, E.; Salzner, U. Tuned Range-Separated Hybrids in Density Functional Theory. Annu. Rev. Phys. Chem. 2010, 61, 85−109. (52) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision B.01; Gaussian, Inc.: Wallingford, CT, 2010. (53) Ito, S.; Nagami, T.; Nakano, M. Density Analysis of Intra- and Intermolecular Vibronic Couplings toward Bath Engineering for Singlet Fission. J. Phys. Chem. Lett. 2015, 6, 4972−4977. (54) Tempelaar, R.; Reichman, D. R. Vibronic Exciton Theory of Singlet Fission. I. Linear Absorption and the Anatomy of the Correlated Triplet Pair State. J. Chem. Phys. 2017, 146, No. 174703. (55) Tempelaar, R.; Reichman, D. R. Vibronic Exciton Theory of Singlet Fission. II. Two-dimensional Spectroscopic Detection of the Correlated Triplet Pair State. J. Chem. Phys. 2017, 146, No. 174704. (56) Tempelaar, R.; Reichman, D. R. Vibronic Exciton Theory of Singlet Fission. III. How Vibronic Coupling and Thermodynamics Promote Rapid Triplet Generation in Pentacene Crystals. J. Chem. Phys. 2018, 148, No. 244701.

15411

DOI: 10.1021/acs.jpcc.9b01713 J. Phys. Chem. C 2019, 123, 15403−15411