Theoretical Study on Third-Order Nonlinear Optical Properties for One-Hole-Doped Diradicaloids.

We investigate the relationships between open-shell character and longitudinal static second hyperpolarizabilities γ for one-hole-doped diradicaloids using the strong-correlated ab initio molecular orbital methods and simple one-dimensional (1D) three-site two-electron (3s-2e) models. As examples of one-hole-doped diradicaloids, we examine H3 +, methyl radical trimer cation ((CH3)3 +), silyl radical trimer cation ((SiH3)3 +), and 1,2,3,5-dithiadizolyl trimer cation (DTDA3 +). For H3 +, the static γ exhibits negative values and shows a monotonic increase in amplitude with an increase in the open-shell character defined by a neighbor-site interaction (y S). On the other hand, it is found for (CH3)3 +, (SiH3)3 +, and DTDA3 + that the static γ value exhibits similar behavior to that for H3 + up to an intermediate y S value, while it takes the negative maximum at a large y S value, followed by a decrease in γ amplitude, and subsequently, γ changes to positive values with a drastic increase for larger y S values. For example, in DTDA3 +, the negative/positive γ values, -69 × 105/700 × 105 au at y S = 0.75/0.87, exhibit significant enhancements in amplitude, 2.4/24 times as large as that (-29 × 105 au) at intermediate y S = 0.59 as is often the case in DTDA2. Using the 1D 3s-2e valence-bond configuration interaction model, these sign inversions and drastic increase in the amplitude of γ are found to originate in the differences in Coulomb interactions between valence electrons, between valence and core electrons, and between valence electrons and nuclei. These results contribute to pave the way for the construction of novel control guidelines for the amplitude and sign of γ for one-hole-doped diradicaloids. Crown

Introduction

Organic nonlinear optical (NLO) molecules have been intensively investigated experimentally and theoretically due to their low cost, high tunability, and fast response.[1] It is known that large second hyperpolarizabilities (γ)—third-order NLO property at the molecular level—of organic molecules can be achieved by optimizing several factors such as π-conjugation size,[2,3] strength of donor/acceptor groups,[4−6] dimensionality of π-conjugations,[7] and amount of hole doping.[8,9] Moreover, Nakano et al. have introduced the open-shell character y (0 (closed shell) ≤ y ≤ 1 (pure open shell)) as an additional factor and have found that symmetric singlet systems with intermediate y (referred to as symmetric singlet diradicaloids) exhibit larger γ than closed-shell or pure open-shell analogues.[10−12] This guideline for third-order NLO molecules based on y has been supported by the experimental observation of gigantic two-photon absorption (TPA) cross section[13,14] and highly efficient third-harmonic generation (THG)[15,16] in several diradicaloids.

To control the third-order NLO properties in diradicaloids, we have also investigated the variations in γ values by changing several physicochemical factors—multiradical character beyond diradicaloids,[17−20] asymmetry electron distribution caused by applying the static electric field[21] as well as by substituting donor/acceptor groups,[22] hole doping,[18,19,23] and the change of spin states.[18−20,24,25] The details of the open-shell NLO design guidelines for symmetric/asymmetric diradicaloids/multiradicaloids as well as of the calculation results of model and realistic open-shell molecular systems are described in our recent review.[26] Among them, hole-doping effects on several diradicaloids have been predicted to show a switching behavior in a sign of γ before and after hole doping. Because negative γ leads to self-defocusing of a light beam, that is, prevention of the damage from the applied strong laser light,[1,28] systems with negative γ are expected to yield durable NLO devices. However, hole-doping effects on diradicaloids have not been revealed in detail.

Therefore, in this study, we investigate the dependence of the static γ on y, referred to as y–γ correlation hereafter, for 1,2,3,5-dithiadizolyl trimer cation (DTDA3+), which is a realistic one-hole-doped diradicaloid, using the strong-correlated ab initio molecular orbital (MO) methods such as the spin-unrestricted coupled cluster singles and doubles (UCCSD) and that with perturbative triples (USSCD(T)) methods.[17,23] As simpler models, we also examine H3+, methyl radical trimer cation ((CH3)3+) and silyl radical trimer cation ((SiH3)3+) models. The origin of the y–γ correlations for these one-hole-doped diradicaloids is revealed within the one-dimensional (1D) three-site two-electron (3s-2e) valence-bond configuration interaction (VBCI) model together with its extended model. The present results contribute to constructing guidelines for controlling the third-order NLO properties in one-hole-doped singlet diradicaloids.

Results and Discussion

Correlation between Pseudo-Diradical Character (yS) and Static Second Hyperpolarizability (γ) for One-Hole-Doped Diradicaloid Models

An open-shell nature for one-hole-doped diradicaloids is evaluated by the diradical character y, which is defined as the occupation number (nLUNO) of the lowest unoccupied natural orbital (LUNO)for whole systems as in the previous study.[23] In addition, we introduce a new open-shell character yS, referred to as “pseudo-diradical character”, originating in a neighbor-site interaction, which is defined by the diradical character (eq ) for the corresponding neutral diradicaloid. For example, the yS values for H3+, (CH3)3+, (SiH3)3+, and DTDA3+ indicate the diradical characters (eq ) for H2, (CH3)2, (SiH3)2, and DTDA2, respectively.

The yS–γ correlations for one-hole-doped diradicaloid models (H3+, (CH3)3+, (SiH3)3+, and DTDA3+, Figure ) as well as for the corresponding neutral diradicaloids (H2, (CH3)2, (SiH3)2, and DTDA2) are shown in Figure . For the neutral diradicaloids, the γ values are positive for any yS and show the maxima (γmax) in the intermediate yS region (γmax = 2.36 × 103 au at yS = 0.56 for H2, γmax = 14 × 103 au at yS = 0.47 for (CH3)2, γmax = 157 × 103 au at yS = 0.46 for (SiH3)2 and γmax = 122 × 103 au at yS = 0.59 for DTDA2). These features in the neutral diradicaloids just follow the y–γ correlation in our previous studies on neutral diradicaloids.[11,12] For the one-hole-doped diradicaloids, γ for H3+ (Figure a) exhibits negative values for any yS and shows a monotonic increase in amplitude with an increase in yS, as shown in a previous study.[23] On the other hand, a realistic system (DTDA3+, Figure d) also exhibits negative γ values and shows an increase in γ amplitude up to a large yS region (yS < 0.75), followed by the decrease in γ amplitude, and subsequently, the change of a sign of the γ value from negative to positive at a larger yS (yS = 0.82). This intriguing behavior is also observed in the other one-hole-doped diradicaloid models ((CH3)3+ (Figure b) and (SiH3)3+ (Figure c)). In (SiH3)3+, although the change of the sign of γ is not found in the region of R < 5.8 Å (yS < 0.91), this feature is expected to appear in the region of R > 5.8 Å, as shown in the case of using larger basis sets (Figure S2). In contrast to the neutral diradicaloids, it is found that the one-hole-doped diradicaloids exhibit the significantly enhanced amplitudes of γ values (negative) at intermediate yS values (though the yS values are larger than those giving γmax for the corresponding neutral diradicaloids) and, except for H3+, give further enhanced amplitudes of γ values (positive) at larger yS values (Table ).

Figure 1

Molecular structures of H3+ (a), CH3 radical (b), (CH3)3+ (c), SiH3 radical (d), (SiH3)3+ (e), DTDA radical (f), and DTDA3+ (g). R denotes the bond distance between neighbor molecules (or atoms) for each one-hole-doped diradicaloid. A longitudinal direction for each model is taken in the y-axis direction.

Figure 2

yS–γ correlations for H2 and H3+ (CISD/aug-cc-pVDZ for yS and γ) (a), for (CH3)2 and (CH3)3+ (PUHF/6-31G* for yS and UCCSD(T)/6-31G* for γ) (b), for (SiH3)2 and (SiH3)3+ (PUHF/6-31G* for yS and UCCSD(T)/6-31G* for γ) (c), and for DTDA2 and DTDA3+ (PUHF/6-31+G for yS and UCCSD(T)/6-31+G for γ) (d). The accuracies of these calculation methods are shown in Tables S1 and S2.

Table 1

γ Valuesa at Several Characteristic yS Values (Shown in Parentheses at Each Data) for One-Hole-Doped Diradicaloids (H3+, (CH3)3+, (SiH3)3+, and DTDA3+)

yS region	H3+ (105 au)	(CH3)3+ (105 au)	(SiH3)3+ (105 au)	DTDA3+ (105 au)
intermediate	–5.45 (0.56)	–21.3 (0.47)	–22.8 (0.46)	–29 (0.59)
large (γ < 0)	–1380 (0.94)	–640 (0.85)	–16 000 (0.89)	–69 (0.75)
large (γ > 0)	–b	2100 (0.92)	–b	700 (0.87)

The number of significant digits in γ values is determined based on the energy convergences and the amplitudes of applied electric fields, which are different from each other (see Tables S1 and S2).

The “–” implies no data corresponding to “large (γ > 0)” due to negative γ values over the whole yS examined for H3+ and for (SiH3)3+, as shown in Figure a,c. Note that in (SiH3)3+, the change of the sign of γ is found to appear in the region of yS > ∼0.9 in the case of using larger basis sets (Figure S2).

To elucidate the difference in yS–γ correlation between the H3+ model and the other models ((CH3)3+, (SiH3)3+ and DTDA3+), we apply the γ density analysis[9] using the third derivative of electron density with respect to F, ρ(3)(), which is referred to as γ densityFor a simple explanation of γ density analysis, we consider a pair of γ densities composed of positive and negative γ densities. The arrow drawn from positive to negative γ density shows the sign of the contribution to γ determined by the relative spatial configuration between the two γ density values. Namely, the sign of the contribution to γ is positive when the direction of the arrow coincides with the positive direction of the coordinate system, while it is negative in the opposite case. The contribution to γ determined by γ density of the two spatial points is more significant when the distance between them is larger. We here conduct the γ density analysis to elucidate the origin of the yS–γ correlations for the one-hole-doped diradicaloids and to discuss the basis set dependence of the yS–γ correlations. The γ density distributions for H3+, (CH3)3+, (SiH3)3+, and DTDA3+ are shown at intermediate yS (Figure a,d,g,j), at large yS with γ < 0 (Figure b,c,e,h,k) and at large yS with γ > 0 (Figure f,i,l). It is found that the γ density distributions for H3+ are localized only around H atoms for any yS, and the γ values always exhibit negative values due to larger negative contributions of γ density on outer H atoms than positive contributions on the inner H atoms. On the other hand, the γ densities for (CH3)3+, (SiH3)3+, and DTDA3+ have characteristic distributions between the inner and outer molecular units, which give the positive contribution to γ, though the γ densities on the outer molecular units give the negative contribution to γ as in H3+. The variation in γ from negative to positive values with an increase in yS for (CH3)3+, (SiH3)3+, and DTDA3+ is found to be described by the variation in the characteristic distributions between the inner and outer molecular units (with positive contribution), where the characteristic distributions grow up with an increase in yS and finally cover over the γ density distribution on the outer molecular units (with negative contribution). The intriguing yS–γ correlations for (CH3)3+ and (SiH3)3+ are expected to be characterized by only one orbital (2p in the C atom for (CH3)3+ and 3p in the Si atom for (SiH3)3+). The difference in yS–γ correlation between H3+ and the others is thus predicted to originate in whether s-type orbital or p-type orbital is used as a characteristic orbital to describe the yS–γ correlation, where the use of p-type orbitals is speculated to cause the sign inversion of γ in the large yS region. From these discussions, we predict that the realistic one-hole-doped diradicaloids such as (CH3)3+, (SiH3)3+, and DTDA3+ exhibit the intriguing yS–γ correlation, as shown in Figure b–d due to the interaction between molecular units with p-type orbitals, for example, the pancake bond.[19,29−31]

Figure 3

γ density distributions for each one-hole-doped diradicaloid (H3+ (R = 2.2 Å (a), 2.8 Å (b), 3.6 Å (c) using the CISD/aug-cc-pVDZ method), (CH3)3+ (R = 2.8 Å (d), 3.8 Å (e), 4.0 Å (f) using the UCCSD/6-31G* method), (SiH3)3+ (R = 4.0 Å (g), 5.6 Å (h), 5.8 Å (i) using the UCCSD/6-31G* method), and DTDA3+ (R = 3.1 Å (j), 3.4 Å (k), 3.6 Å (l) using the UCCSD/6-31+G method)) at intermediate yS ((a) (yS = 0.56), (d) (yS = 0.47), (g) (yS = 0.46), (j) (yS = 0.60)), at large yS with γ < 0 ((b) (yS = 0.81), (c) (yS = 0.94), (e) (yS = 0.85), (h) (yS = 0.89), (k) (yS = 0.75)), and at large yS with γ > 0 ((f) (yS = 0.89), (i) (yS = 0.91), (l) (yS = 0.82)). The γ values are also shown for each one-hole-doped diradicaloid. The yellow and blue surfaces represent the positive and negative γ density distributions with ±104 au (a, d, j, k, l), ±2 × 104 au (g), ±105 au (b, e, f), ±106 au (h, i), and ±107 au (c), respectively.

One-Dimensional (1D) Three-Site Two-Electron (3s-2e) Valence-Bond Configuration Interaction (VBCI) Models

One-Dimensional 3s-2e VBCI Scheme

The valence-bond configuration interaction (VBCI) models have been employed to investigate hyperpolarizabilities for charge-transfer closed-shell organic molecules,[5] and symmetric[12] and asymmetric neutral diradicaloids.[21,22] However, these VBCI models have rarely been extended to multisite (multiradical) models involving more than two sites. On the other hand, the giant third-order NLO responses in inorganic materials composed of 1D charge-transfer Mott insulators have been reported[33] and analyzed using the extended Hubbard models such as single-band[34] and two-band Hubbard models,[35] where experimental data have been used as the parameters in these models. Note here that their models have never treated the open-shell (multiradical) singlet systems with a wide range of open-shell characters and have never been used to clarify the relationships between open-shell character, structure, and γ, for example, the bond distance dependence of γ through a variation in open-shell character. Thus, we construct the VBCI model for the one-hole-doped diradicaloids, which is referred to as 1D 3-site 2-electron (3s-2e) VBCI model, to clarify the origin of the intriguing yS–γ correlations with changing bond distances. We start from formulating the 1D 3s-2e VBCI model.

As the simplest example of one-hole-doped diradicaloids, we consider a cationic three-site model (A1–A2–A3)+ with two electrons in three orbitals (a1, a2, and a3 on A1, A2, and A3, respectively, referred to as VB one-electron orbital (VB-OEO) for each orbital), in which three sites (A1, A2, and A3) are the same as each other, such as H3+, (CH3)3+, (SiH3)3+, and DTDA3+ (Figure ). Each VB-OEO is localized around each site, corresponding to singly occupied molecular orbital (SOMO) of monomers for these systems. Note that these orbitals are normalized but not orthogonalized to each other so thatFor MS = 0 (singlet and triplet), using VB-OEOs, there are six covalent, |a1a̅2⟩ (≡|a1a̅2core⟩), |a2a̅1⟩ (≡|a2a̅1core⟩), |a1a̅3⟩ (≡|a1a̅3core⟩), |a3a̅1⟩ (≡|a3a̅1core⟩), |a2a̅3⟩ (≡|a2a̅3core⟩), and |a3a̅2⟩ (≡|a3a̅2core⟩), and three ionic, |a1a̅1⟩ (≡|a1a̅1core⟩), |a2a̅2⟩ (≡|a2a̅2core⟩), and |a3a̅3⟩ (≡|a3a̅3core⟩), determinants, where “core” means the closed-shell core orbitals orthogonal to other core orbitals and to VB-OEOs (⟨c|c′⟩ = 0,⟨a|c′⟩ = 0(c,c′ ∈ core,i = 1,2,3)), and upper bar (nonbar) indicates the β (α) spin. From the definition of ab initio N-electron Hamiltonianwhere N indicates the number of electrons in VB-OEOs and core orbitals. The CI matrix in the VB-OEO representation {|a1a̅2⟩, |a2a̅1⟩, |a1a̅3⟩, |a3a̅1⟩, |a2a̅3⟩, |a3a̅2⟩, |a1a̅1⟩, |a2a̅2⟩, |a3a̅3⟩} takes the formwhere a form after the equal sign expresses only the upper triangular elements. The Coulomb interaction between neighbor sites (U (≡⟨a1a2|a1a2⟩)) subtracted from the energy of the |a1a̅2⟩ configuration (⟨a1a̅2|Ĥ|a1a̅2⟩) is taken as the energy origin. Here, the notations in this CI matrix are defined bywhere, in eq where the off-diagonal elements in eq ignore the terms of squared overlap integrals (S2 ∼ 0) since they are significantly small in the region (yS > 0.5) treated in the first approximation explained in the next paragraph. The derivations of all of the matrix elements in eq are shown in the Supporting Information. To compare the yS dependence of γ values obtained from the diagonalization of eq including 15 variables with that obtained from ab initio calculations, we reduce these variables in eq to only one variable. However, the analytical relations of two-electron integrals in eq cannot be obtained except for that composed of only 1s-type Gaussian-type orbitals as VB-OEOs. Thus, we first formulate the matrix elements in the case of H3+, which are composed of only 1s orbital and an empty of core orbitals, and then construct the 1D 3s-2e VBCI model for the other one-hole-doped diradicaloids ((CH3)3+, (SiH3)3+ and DTDA3+) by including the rest of contributions such as the expansion effects of electron distributions in core orbitals and VB-OEOs.

Figure 4

Schematic diagram of the 1D 3s-2e VBCI model (with site labels (A1, A2, and A3) and VB one-electron orbitals (VB-OEOs, a1, a2, and a3)) composed of covalent and ionic configurations. A longitudinal direction in this model is taken in the y-axis direction. The up (red) and down (blue) arrows indicate valence electrons with up and down spins, respectively. This model is composed of nine configurations with MS = 0, which are six covalent configurations (in yellow panel) and three ionic configurations (in green panel), and provides six singlet and three triplet states by solving the eigenvalue equation with the 1D 3s-2e VBCI model Hamiltonian. Here, the pseudo-diradical character yS is predicted to change with the intersite distance R.

One- and two-electron integrals between 1s-type Gaussian-type orbitals with a common orbital exponent α are expressed by[36]where and indicate a position vector of a center of orbital a and that of a nucleus, respectively; Z and r indicate a nuclear charge at site A, +1 for H3+ model, and a distance between an electron and the nucleus at site A, respectively; P (Q) indicates a position vector at the middle of and ( and ). Because we are interested in the region of large yS, namely, weak neighbor-site interactions, the overlap integrals between VB-OEOs on different sites are predicted to become significantly small values. Therefore, the terms proportional to S13 or squared overlap integral are assumed to be negligible after converting each integral into the expression using the overlap integrals, leading toEquations –20 are the first approximations of all.

Second, we adopt the following relationswhere R indicates | – | and Z is equal to +1 for H3+. These relations imply that an electron in a VB-OEO (a) is considered to be completely localized at site A, as shown in the Supporting Information, and can be sufficiently satisfied in the case of a large R, that is, a large yS.

Third, we set the common values (U and t12) as the on-site Coulomb interactions (U) and the neighbor-site transfer integrals (t12() with respect to any i, respectivelyThe former relation (eq ) is exactly equal to each other when α is the same for all of the VB-OEOs, whereas the latter one (eq ) corresponds to the definition of transfer integral for 2s-2e VCI model[12] and is not exact but can provide good approximations to the excitation energies and transition moments for H3+ at large R (yS), as shown in Figure S6.

Fourth, we express the neighbor-site transfer integral, t12, as a function of the neighbor-site Coulomb interaction, U. From eqs –15, t12 can be considered to be roughly proportional to neighbor-site overlap integral (S12; see also eq ), that isBy substituting eq into eq , t12 is approximately expressed as a function of Uwhere A and B are constants. Although the constant B in eq can be estimated to α/2 from eq , the proportional relation (eq ) is an approximation so that the B slightly deviates from α/2. Therefore, the constant B should be determined by fitting t12 to U for H3+.

By applying the above four approximations (eqs –20, 21 and 22, 23 and 24, and 26) to eqs , 11 variables in the VBCI matrix are reduced to two variables, U and U. Here, we define two dimensionless quantitiesThe VBCI eigenvalue equation is expressed bywhere corresponds to the Hamiltonian matrix (eq ). Dividing by on-site Coulomb interaction (U), this equation is rewritten into the dimensionless formwhere DL (DL ≡ /U) and DL (DL ≡ /U) indicate dimensionless Hamiltonian and eigen energy matrices, respectively. Using eqs and 28, the explicit expression of DL after applying the four approximations is given byNote that h = 1/2U, as shown in the Supporting Information. This dimensionless Hamiltonian DL includes one variable r due to the relation between r and r defined in the fourth approximation (eq ). Therefore, the results obtained from solving the eigenvalue equation (eq ) with the dimensionless Hamiltonian (eq ) are expected to correspond one to one to the results of the H3+ model. Here, the model described by the dimensionless Hamiltonian (eq ) is called the 1D 3s-2e VBCI model.

To describe the open-shell character between neighbor sites, we introduce the pseudo-diradical character (yS) defined in two-site diradical systemsthe derivation of which is shown in the Supporting Information. Note that “U” involved in r and r in eq indicates the on-site Coulomb interaction, which has a different definition from that in the previous study.[12,21,22]

Next, we obtain the concrete relation formula between r and r for H3+. The parameters for H3+ should be obtained by CISD/STO-1G calculation using GAMESS program package[37] due to one 1s-type Gaussian-type orbital per site in the 1D 3s-2e VBCI model. It is found that the yS–γ correlation for H3+ using the STO-1G basis set qualitatively reproduces that obtained using larger basis sets, as shown in Figures a and 5. The on-site Coulomb interaction (U) is a constant so that the relationship between r and r is expressed byHere, A′ and B′ are defined by A′ ≡ U–1A and B′ ≡ U–2B, respectively, both of which are constants. By fitting the relationship between r and r using eq for H3+ with the STO-1G basis set (Figure S5), A′ and B′ are obtained as shown in Table . U = 0.728 is obtained from the calculation of the two-electron integral for STO-1G orbital of H atom, which is a in eq (i = j) so that B = 0.18. In GAMESS program package, an orbital exponent of 1s orbital for H atom (α) is equal to 0.4166 so that the approximate parameter B, B ∼ α/2 as mentioned above (eqs and 26), is equal to 0.2083. It is found that this approximate value (0.2083) obtained from the orbital exponent is close to a value (0.18) estimated by eq .

Figure 5

yS–γDL correlations for the 1D 3s-2e VBCI and H3+ models. The γDL values for the H3+ model are calculated using eq , in which U = 0.728 and R in atomic units are adopted and γ values are calculated by the perturbative expression (eq ) using the CISD/STO-1G results.

Table 2

Fitting Parameters in Equation for H3+ Using the STO-1G Basis Set

	A′	B′
values	1.50	0.34

According to the perturbation theory, the diagonal static γ along the y-axis, which corresponds to a direction of aligned three sites, for the centrosymmetric systems are expressed bywhere μ and E indicate the transition moment between the ath and bth excited states (a = 0 for the ground state) and excitation energy of the ath excited state, respectively. Note that γ is expressed in the B-convention. In the 1D 3s-2e VBCI model, the dimensionless γ (γDL) is used instead of γ given in eq as in the previous study on the symmetric and asymmetric 2s-2e V(B)CI models,[12,21,22] and can be defined bywhere R12 indicates the distance between the neighbor sites; μDL  and EDL indicate dimensionless transition moment between eigenstates a and b, and dimensionless excitation energy of excited state a, respectively, defined byThe transition-moment matrix, μ, in VB-OEO representation is defined bywhere r1 and r2 indicate the positions of two electrons in VB-OEOs. Here, we assume that the coordinate origin is set to site A2 (see Figure ), and then μ13,13 = μ22,22 = 0, so that the diagonal elements of the transition-moment matrix (eq ) are expressed byof which the derivation is shown in the Supporting Information. In addition, we adopt the approximation that the off-diagonal elements are ignored. In the 1D 3s-2e VBCI model, this approximation is found to be sufficient to reproduce the transition moment between two states for the H3+ model quantitatively, as shown in Figure S6b in the region from intermediate yS to pure open shell (yS > 0.5). The eigenvectors and eigenvalues obtained by diagonalization of the 1D 3s-2e VBCI matrix (eq ) give the dimensionless transition moments and the excitation energies, respectively, and, using these properties, the γDL (eq ) value can be calculated.

Comparison between 1D 3s-2e VBCI Model and H3+ Model

In this subsection, we compare the yS–γDL correlation for the 1D 3s-2e VBCI model with that for the H3+ model to verify the validity of the 1D 3s-2e VBCI model. The yS–γDL correlations for the 1D 3s-2e VBCI and H3+ models are shown in Figure . Note that γDL defined in eq for the H3+ model is used, in which U = 0.728 and R in atomic units are adopted and γ values are calculated by the perturbative expression (eq ) using the CISD/STO-1G results. The yS–γDL correlation for the 1D 3s-2e VBCI model is found to coincide well with that for the H3+ model, namely, the four approximations (eqs –20, 21 and 22, 23 and 24, and 26) and a treatment of ignoring the off-diagonal elements in eq are reasonable in the region from intermediate yS (yS ∼ 0.5) to pure open shell (yS ∼ 1). This fact is also clear from the relationship between yS and the dimensionless excitation energies (Figure S6a), and between yS and the dimensionless transition moments (Figure S6b) for the 1D 3s-2e VBCI and H3+ models.

Radical Site Character (χ) Dependences of Correlation between Pseudo-Diradical Character (yS) and Dimensionless Second Hyperpolarizabilities (γDL)

To clarify the origin of the intriguing yS–γ correlations for (CH3)3+, (SiH3)3+ and DTDA3+ models (Figure ), we modify the 1D 3s-2e VBCI model Hamiltonian including the rest of contributions such as the expansion effects of electron distributions in the core orbitals and VB-OEOs. The (1,1), (2,2), (5,5), and (6,6) elements in eq are the same due to the centrosymmetric system, but these and (3,3) or these and (4,4) are not exactly equal to each other. The equivalence for these six elements in eq arises from considering each VB-OEO (a) as completely localized at each site (A) (see the Supporting Information). In the H3+ model, this treatment is found to be sufficient to quantitatively reproduce yS–γ correlation for ys > 0.5 because the delocalization in electrons of 1s-type orbital used as the VB-OEO can be negligible. In other models ((CH3)3+, (SiH3)3+, and DTDA3+), however, the VB-OEOs become p-type orbital with electron distribution significantly expanding toward the aligned three sites so that each VB-OEO cannot be considered to be completely localized at each site.

To include the expansion effects of p-type and of the core orbitals, we introduce the “pseudo-classical Coulomb interaction model”, where the VB-OEOs (a, i = 1–3) and the core orbitals with similar expansion to the VB-OEOs (referred to as sub-VB-OEOs) are assumed to have certain expansions (Figure ). For (CH3)3+, the VB-OEO and sub-VB-OEOs on each site, which correspond to 2p orbital, and 2p, 2p, and 2s orbitals, respectively, have certain expansions, whereas the expansions of 1s orbitals are ignored. In this case, the expansions of 2s and 2p orbitals are estimated from the “radical site characters of 2s and 2p orbitals (χO, O = s, p)” defined by the ratio of the effective size (Rs and Rp for 2s and 2p orbitals, respectively) to the bond distance (R12)where χO = 0 corresponds to the 1D 3s-2e VBCI model ignoring the expansions of VB-OEOs and sub-VB-OEOs. In the pseudo-classical Coulomb interaction model, it is assumed that the electron is distributed at the position of a distance of RO (O = s, p) from the sites with an equal probability so that, for (CH3)3+, the terms in diagonal elements of the Hamiltonian matrix in eq are expressed bywhere, in eq These derivations of eqs –46 are shown in the Supporting Information. Note that these functions, fnuc,χ, fp,χ, and fs,(χ, depend on R12 through χO (O = s, p). Using these relations, the dimensionless Hamiltonian (eq ) for (CH3)3+ are rewritten in the following form including the expansions of VB-OEOs and sub-VB-OEOswhere χ indicates a net total radical site character, defined bywhich is a dimensionless parameter as a function of χs and χp and takes positive values as shown in Figure in the next paragraph. Hence, χr can be understood as a relative destabilization of valence electrons (one-electron dimensionless energy) on site A2 with respect to those on sites A1 and A3 (Figure a) or as a relative stabilization of centrosymmetric covalent valence-electron configurations with respect to those of asymmetric covalent valence-electron configurations (Figure b). These relative energy variations caused by considering the distribution expansions of valence electrons and core electrons are understood by the enhancement of the Coulomb repulsion from electrons localized on the neighbor sites, stabilizing the symmetric covalent configuration states compared to the asymmetric covalent ones. Although the (8,8) element in eq , the energy level of ionic configuration |a2a̅2⟩, should be modified, the ionic configurations are predicted to give a negligible contribution to γDL in the region from intermediate yS to pure open shell (yS > 0.5) due to high-lying energy levels for the ionic configurations. Thus, we ignore the radical site character for the (8,8) element. Also, for (SiH3)3+, the dimensionless Hamiltonian matrix obtained in the similar approximations is the same form as that for (CH3)3+, while that for DTDA3+ is different from those for (CH3)3+ and (SiH3)3+ due to the difference in the number of the sub-VB-OEOs, which include p-type sub-VB-OEOs with expansions toward the aligned three sites. However, from the intriguing yS–γ correlation for (CH3)3+ and (SiH3)3+ similar features to that for DTDA3+, this difference for DTDA3+ can be considered to be emerged as differences in χ, and that in the fitting parameters in eq , so that the same form of the dimensionless Hamiltonian matrix for DTDA3+ can be obtained. Because the intriguing yS–γ correlation is predicted to originate in resolving the energy degeneracy of six covalent configurations, the relationship between χr and the off-diagonal element (r) is predicted to play an important role of the intriguing yS–γ correlation. Namely, if the fitting function for r and r has the different parameters (A′ and B′) in eq , the intriguing yS–γ correlation can be reproduced qualitatively by adjusting χ appropriately. Therefore, we use the fitting function with parameters for H3+ in Table , and the missing effects in the 1D 3s-2e VBCI model corresponding to the H3+ model are included by adjusting the χ value. In addition, the intriguing yS–γ correlation for DTDA3+ is also observed in the case of (CH3)3+ and (SiH3)3+, as shown in Figure , so that we elucidate the intriguing yS–γ correlation using the extended dimensionless Hamiltonian matrix with two-electron two-center bonding between neighbor sites described by the pseudo-classical Coulomb interaction model as in Figure (eq ).

Figure 6

Coulomb interactions between two VB-OEOs (2p orbitals) (a), between VB-OEO and nucleus (b), between VB-OEO and s-type sub-VB-OEO (2s orbital) (c), and between VB-OEO and p-type sub-VB-OEO (2p (d) and 2p (e)) in the pseudo-classical Coulomb interaction model. The electron in VB-OEO or p-type sub-VB-OEO is considered to be localized at the position of a distance of Rp away from the site A, that is, the electron in VB-OEO is localized at these two positions each with a probability of 50%. On the other hand, the electron in s-type sub-VB-OEO is considered to be uniformly distributed on the sphere of the distance of Rs from site A.

Figure 8

(χs, χp) dependences of χ (eq ). The ranges of both χs and χp are [0, 0.4], and the color range of χ goes from 0.0 (violet) to 0.2 (red).

Figure 7

Relative one-electron (a) and two-electron dimensionless energies (b) before and after considering the expansions of VB-OEOs and sub-VB-OEOs. The broken lines indicate these energies before considering the expansion effects in the 1D 3s-2e VBCI model. “Symmetric Covalent” and “Asymmetric Covalent” indicate the symmetric and asymmetric covalent valence-electron configurations, respectively.

The (χs, χp) dependences of χ are shown in Figure . Note that both χs and χp take only positive values, and the expansion of 2p orbital is larger than that of 2s orbital (χs < χp). It is found that χ exhibits a positive value for any (χs, χp) and increases with increases in (χs, χp). This fact provides the lower energy level for the symmetric covalent valence-electron configurations (|a1a̅3⟩ and |a3a̅1⟩) than that for the asymmetric covalent valence-electron configurations (|a1a̅2⟩, |a2a̅1⟩, |a2a̅3⟩, and |a3a̅2⟩) as shown in Figure b, and then the symmetric configurations become dominant in the ground state at large yS, while the first and second excited states are mainly composed of the asymmetric configurations. Such electronic structures for these eigenstates are similar to that for the neutral symmetric diradicaloids (2s-2e systems) with positive γ values for any yS, in which the symmetric and asymmetric configurations correspond to the covalent and ionic configurations in the previous study.[38]

The yS–γDL correlation for each χ value is shown in Figure . Although χ decreases with an increase in yS because χO (O = s, p) decreases with an increase in R12 from eq due to constant RO (O = s, p) for any R12, and yS increases with an increase in R12 from eqs and 32, we use a constant value of χ for any yS. In spite of this treatment, surprisingly, the yS–γDL correlation for this model is found to qualitatively reproduce the features of the intriguing yS–γ correlations for (CH3)3+, (SiH3)3+, and DTDA3+ (see Figures and 9). As seen from Figure , yS giving the negative extremum of γDL slightly decreases with an increase in χ, the amplitude of the negative extremum of γDL is significantly suppressed with an increase in χ, and yS changing the sign (from negative to positive) of γDL decreases with an increase in χ. As mentioned above, χ involves the effect of the difference in the fitting function between H3+ and the others ((CH3)3+, (SiH3)3+, and DTDA3+). If the result of the yS–γ correlation in Figure completely matches that of this model, χ values are roughly estimated at 0.04 for (CH3)3+ and (SiH3)3+, and 0.08 for DTDA3+ by judging from yS giving the negative extremum of γ. These differences in χ are predicted to originate from the difference in whether including the p-type orbitals as sub-VB-OEOs. As seen from the first three terms for sub-VB-OEOs in eq , the expansions of sub-VB-OEOs contribute positively to the χ value. For (CH3)3+ and (SiH3)3+, sub-VB-OEOs are composed of only s-type, p-type, and p-type orbitals, while sub-VB-OEOs for DTDA3+ additionally include p-type orbitals so that the χ value for DTDA3+ becomes larger than those for (CH3)3+ and (SiH3)3+ due to larger expansions of p-type orbitals than that of the other type orbitals.

Figure 9

yS–γDL correlations for χ = 0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, and 0.20. yS–γDL correlations in the range [0, 5000] of γDL are shown in the inset.

On the other hand, the intriguing yS–γ correlation for (CH3)3+, (SiH3)3+, and DTDA3+ originate in resolving the energy degeneracy of six covalent configurations due to χr as in Figure so that r is one of the key factors and depends on the expansions of s-type and p-type orbitals as well as conjugated size of the molecule in the previous study for organic metals.[39] DTDA radical has larger π-conjugation than (CH3) and (SiH3) radicals so that r for DTDA is predicted to be smaller than those for (CH3)3+ and (SiH3)3+ due to the reduction of the neighbor-site Coulomb interaction. In spite of this fact, χr for DTDA3+ is found to be larger than those for (CH3)3+ and (SiH3)3+. Namely, this is predicted to be caused by the fact that the increasing effect of χ exceeds the decreasing effect of r in DTDA3+. In case of larger radical molecules than DTDA, however, the r value can be more reduced so that the decreasing effect of r may exceed the increasing effect of χ in a specific molecular size, resulting in the reduction of χr. To verify this prediction, more detailed investigation on the molecular size dependence of yS–γ correlation will be needed in the future.

Note that these results are obtained under the four approximations for H3+ (eqs –20, 21 and 22, 23 and 24, and 26) based on the fitting parameters in eq for H3+ using STO-1G basis set (Table ), as well as based on modifying the diagonal elements in eq by the pseudo-classical Coulomb interaction model. In addition, configurations considered in the extended 1D 3s-2e VBCI model are restricted to only nine configurations. Despite applying these approximate treatments, the obtained yS–γDL correlation is found to qualitatively reproduce the yS–γ correlation obtained by the strong-correlated ab initio MO methods so that the present VBCI model is expected to be one of the promising methods to analyze yS–γ correlation for one-hole-doped diradicaloids as well as general Ms-Ne systems, where M and N indicate the number of sites and electrons, respectively. In particular, because 1D one-hole-doped multiracaloids (Ms-(M – 1)e systems) are expected to exhibit the intriguing yS–γ correlations similar to the one-hole-doped diradicaloid, elucidating that the yS–γ correlations for Ms-(M – n)e (n ≥ 1) systems will be a challenging research topic in the next stage.

Conclusions

We investigate the relationship between the open-shell character (yS) and the longitudinal static second hyperpolarizabilities (γ)—yS–γ correlation—for one-hole-doped diradicaloids, H3+, (CH3)3+, (SiH3)3+, and DTDA3+, and the corresponding neutral diradicaloids, H2, (CH3)2, (SiH3)2, and DTDA2, by the UCCSD and UCCSD(T) methods. For the neutral diradicaloids, γ exhibits positive values for any yS and a maximum value at the intermediate yS for all of the systems as shown in the previous studies.[9−11] On the other hand, for H3+, γ exhibits negative values for any yS and increases in amplitude with an increase in yS, as shown in the previous study,[23] while the yS–γ correlations for (CH3)3+, (SiH3)3+, and DTDA3+ are found to be different from that for H3+: γ (with negative sign) increases in amplitude with an increase in yS up to an intermediate yS, and it takes the negative maximum at intermediate yS, followed by a decrease of γ amplitude, and subsequently, the γ value changes to positive with a drastic increase at a larger yS. Moreover, it is found that all of the one-hole-doped diradicaloids exhibit remarkably enhanced amplitudes of γ (with negative sign) at intermediate yS compared to the neutral diradicaloids. For example, the γ amplitude for DTDA3+ (at yS = 0.75) exhibits 2.9 times as large as that for s-indaceno[1,2,3-cd;5,6,7-c′d′]diphenalene (IDPL),[40] which is a typical neutral diradicaloid with a similar yS value well known as one of the largest TPA cross sections in the similar-sized organic molecules.[13] From γ density analysis, it is predicted that the difference in the yS–γ correlations between H3+ and the others originates from the difference of types of orbitals involved in the interaction between the neighbor sites. In addition, we analyze the difference in yS–γ correlations using the one-dimensional (1D) three-site two-electron (3s-2e) valence-bond configuration interaction (VBCI) model. We formulate the 1D 3s-2e VBCI model and demonstrate that this model reproduces quantitatively the yS–γDL correlation for H3+. Furthermore, we extend the 1D 3s-2e VBCI model by including “radical site characters”, which represent the expansions of valence-bond one-electron orbital (VB-OEO), corresponding to the 2p orbital on each site for (CH3)3+, and the core orbitals similar to VB-OEOs (sub-VB-OEO), corresponding to 2s, 2p, and 2p orbitals on each site for (CH3)3+, and show that the yS–γDL correlation in this extended model provides the behaviors similar to those in (CH3)3+, (SiH3)3+, and DTDA3+. It is turned out from the analysis of this model that the difference in yS–γ correlations between H3+ and the others originates in the difference in the expansions of VB-OEOs and sub-VB-OEOs. These results demonstrate that the one-hole-doped diradicaloids become one of the promising candidates for building blocks of novel third-order NLO materials, the γ amplitudes and their signs of which can be sensitively controlled by subtle chemical modifications. Indeed, the hole-doped aggregates composed of thiazyl radicals have been already synthesized[41] by doping counteranion. Although these aggregates form infinite columnar stacks in the crystalline state, it is expected that the electric structures similar to one-hole-doped diradicaloids are realized using organic-pillared coordination cages,[42] for example. Additionally, the present models are also predicted to offer the new promising tools for clarifying the yS–γ correlation for M-site N-electron (Ms-Ne) systems, such as the realistic open-shell molecular aggregates composed of cyclic thiazyl radicals,[20,43] pentadienyl radicals,[30] phenalenyl radicals,[19,31] and olympicene radicals.[32]

Recently, several groups have discussed the nature of the third-harmonic scattering (THS), which is an incoherent process and originates in γ(−3ω;ω,ω,ω) at the molecular scale, in an isotropic medium like a liquid and a solute in a solution due to its applicability in technique to monitor molecules, ionic species, and so on.[44−46] The γ values related to THS can be divided into three components, isotropic, quadrupolar, hexadecapolar components, which have been observed as the difference in the THS response.[45] From the definition, for 1D one-hole-doped DTDA clusters, the isotropic and quadrupolar components are predicted to be primary components. Because the 2D and 3D hole-doped multiradicaloids are expected to exhibit intriguing dependences of these components on yS, an extension to multidimensional multiradicaloid structures will be also one of the interesting topics.

Computational Methods

The geometry optimization was performed using the UCCSD(T)/aug-cc-pVDZ method for CH3 and SiH3 radicals under the constraint of D3 symmetry and using the UMP2/6-311+G* method[20] for DTDA radical under the constraint of C2 symmetry. DTDA and CH3 radicals were confirmed to be in the stable local minima by the vibrational analysis, while a SiH3 radical was confirmed to be in the saddle point under the constraint of D3 symmetry. The fact that the SiH3 radical exhibits the nonplanar stable local minimum is understood by a difficulty of hybridization of the 3s and 3p orbitals.[47−49] However, since we are interested in the π–π stacking interactions between neighbor radical sites with p-type orbitals, the planar SiH3 radical model in the saddle point under the constraint of D3 symmetry was adopted in this study. The 1D face-to-face stacking trimers are composed of these optimized monomers with a stacking distance R. Also, the geometry of H3+ was built by aligning three hydrogen atoms with a bond distance R. The distance R is varied from 1.2 to 4.0 Å for H3+, from 2.0 to 4.4 Å for (CH3)3+, from 3.0 to 5.8 Å for (SiH3)3+, and from 2.8 to 3.8 Å for DTDA3+. These distance variations correspond to the wide range of yS variations: yS = 0.08–0.97 for H3+, 0.09–0.94 for (CH3)3+, 0.02–0.91 for (SiH3)3+, and 0.36–0.87 for DTDA3+.

The yS values were calculated using the configuration interaction singles and doubles (CISD)/aug-cc-pVDZ method for H3+, the PUHF/6-31G* method for (CH3)3+ and (SiH3)3+, and the PUHF/6-31+G method for DTDA3+, in which PUHF is an abbreviation of the spin-projected (P) unrestricted Hartree–Fock (UHF) and yS at the PUHF level is obtained from the occupation number (nHONO) of the highest occupied natural orbital (HONO) and nLUNO at the UHF level by[50,51]The diradical character at the PUHF level is known to well reproduce those using highly correlated CI methods.[51]

Because we focus on the effects of π–π stacking interaction on the electronic structures and properties, only longitudinal components of static γ (γ) in the stacking direction (y) (simply called γ hereafter unless specified otherwise) were evaluated. The γ values were calculated using the finite field (FF) method[52] with B-convention,[53] that is, using the fourth-order numerical differential formula of energy E with respect to applied electric fields Fwhere E(F) indicates the total energy in the presence of F. Each energy in the FF method was calculated using the CISD/aug-cc-pVDZ method for H3+, the UCCSD(T)/6-31G* method for (CH3)3+ and (SiH3)3+, and the UCCSD(T)/6-31+G method for DTDA3+. We also calculated the static γ values for the corresponding dimers (H2, (CH3)2, (SiH3)2, and DTDA2) using the same method as in these trimer cations. The numerical accuracies for the γ calculated by the FF method are shown in Tables S1 and S2. Also, the γ values using these basis sets were confirmed to qualitatively reproduce those using larger basis sets, respectively (see Figures S1 and S2).

γ density is also calculated using the third-order numerical differential formula of electron density (ρ(,F)) with respect to FThe γ value is expressed byHere, r denotes the longitudinal (y) component of the electronic coordinate . Due to the unavailability of calculating electronic density at the UCCSD(T) level of theory by Gaussian09, we perform the γ density analysis at the UCCSD level of theory for (CH3)3+, (SiH3)3+, and DTDA3+. This method is expected to give a reliable feature of γ density distribution since the UCCSD method provides the qualitatively similar behavior of yS–γ correlation (Figure S3). All of these calculations were performed by the Gaussian09 program package.[54]