How to get deeper insights into the optical properties of lanthanide systems: a computational protocol from ligand to complexes

2.1 Ground state geometries for L, EuL₃E₂ and EuL₃T₂

Positional disorder in X-ray crystal structure (CCDC 1539534) suggests that L occurs as two rotamers (Fig. S1; hereafter, A and B) whose percentages amount to 85 ± 0.5 % (A) and 15 ± 0.5 % (B) [35]. In A, both thiophene rings “point” in the same directions, with two S atoms on the same side of the diketone oxygen atoms, while in B the S atoms point at opposite sides. Geometrical parameters of both A and B have been optimized by adopting either the PBE or the B3LYP functionals and, independently of the adopted exchange-correlation (XC) functional, A resulted more stable than B (0.39 kcal/mol and 0.50 kcal/mol in PBE and B3LYP calculations, respectively). All computational details for calculations are reported in Supplementary Information. If a Boltzmann population is assumed for the two rotamers and the energy differences are taken into account, the theoretical distribution at room temperature (RT, 298.15 K) can be obtained by:

where P i is the population of the i-th rotamer, ε i is the energy of the i-th rotamer, k B is the Boltzmann constant and T is the temperature. B3LYP and PBE distributions (69.9 %/30.1 % and 65.9 %/34.1 %, respectively) are very similar and in satisfactory agreement with crystallographic data (optimised structures perfectly match the crystallographic geometries (Fig. S2)), with tiny better accordance for B3YLP. The small difference between theoretical and experimental distributions may be due to crystal packing effects, which cannot be modelled in single molecule calculations. Before going on, it deserves to be emphasized that a third rotamer C with both thiophene rings pointing in the same direction but with the S atoms on the opposite side of the diketone oxygen atoms (Fig. S1) could be present. Even though C has not been experimentally revealed, it has been included in our analyses. Indeed, its absence in the solid state does not necessarily rule out its possible presence in the solution. As such, C is less stable than A by 1.83 kcal/mol and the B3LYP Boltzmann populations of A, B and C correspond to 68 %, 29 % and 3 %, respectively. To verify whether different forms may undergo interconversion, a linear transit calculation has been carried out. At RT, only a negligible fraction (corresponding to 10⁻⁹ of the molecules) would possess sufficiently high energy to overcome the 8 kcal/mol barrier; thus, possible interconversion phenomena have been neglected (details on this matter are reported in the Supplementary Information (Fig. S3)). Even though solvent effects have been considered for A, it has to be underlined that the corresponding B3LYP and PBE optimised structures are perfectly superimposable to the gas-phase ones, thus confirming that the geometry optimisation of small organic molecules experimentally found either in the crystal phase or dissolved in an apolar solvent can be confidently carried out by exploiting gas phase calculations [36, 37].

X-ray data for EuL ₃ E ₂ (CCDC 1539532) and EuL ₃ T ₂ (CCDC 1539535) are reported in ref [19, 35] and corresponding geometries have been used as a starting point for the geometry optimisation calculations. To simplify the comparison between theoretical and experimental structures, the following key quantities have been monitored throughout the numerical experiments: (i) internuclear distances between Eu³⁺ and the coordinated O atoms of L; (ii) the angle formed by the O atoms of L and Eu; (iii) the angle formed by Eu and the central C atoms of two distinct L (see Fig. S4). In fact (i) and (ii) provide information about the Eu³⁺ coordination environment, while (iii) supply details about the L relative orientation. The comparison of the EuL ₃ E ₂ and EuL ₃ T ₂ optimized geometrical parameters with the experimental outcomes (see Tables S1 and S2, respectively) reveals the very good agreement between experiment and theory (Fig. S5), with the Eu-O distances slightly overestimated by ∼2 % on average (see Table S3). Again, tiny deviations may be tentatively attributed to crystal packing phenomena. Before moving to the analysis of spectroscopic data, let us emphasise that B3LYP and PBE geometrical parameters are very similar (Fig. S6), with the former performing slightly better than the latter but at a significantly higher computational cost: the EuL ₃ E ₂ B3LYP geometry optimisation took 17 times more CPU time than the PBE one with the same starting geometry. As the differences between B3LYP and PBE optimized geometries are almost negligible, EuL ₃ T ₂ calculations have been carried out by exploiting only the PBE functional (the PBE optimised geometry of EuL ₃ T ₂ is compared with its X-ray structure in Fig. S7). As a whole, DFT calculations provide a reliable prediction of the ground state geometry of the Eu³⁺ complexes and the optimised structures are in very good agreement with X-ray data.

2.2 L, EuL₃E₂ and EuL₃T₂ absorption spectra

The L absorption spectrum has a single broad band in the near UV region (around 375 nm) and two evident shoulders at 395 and 360 nm (Fig. 2, solid line).

Fig. 2:

Normalized experimental absorption spectrum of L (solid black line) recorded in solution (toluene) superimposed to the different theoretical vertical excitation energies on the rotamer A estimated by adopting the B3LYP optimised ground state geometry but diverse XC functionals (SAOP, LB94, PBE0, B3LYP) (a). 3D plots of the MOs mainly contribute to the electronic transition (b).

TD-DFT numerical experiments have been carried out by exploiting diverse XC functionals (SAOP, LB94, PBE0 and B3LYP) and by using the B3LYP optimized ground state geometry. Regardless of the adopted XC functional, the near UV region includes only one intense transition (Fig. 2a and Table S4) involving the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbitals (LUMO), both of them delocalized over the whole molecule (Fig. 2b). PBE0 and B3LYP XC functionals overestimate the HOMO–LUMO ΔE, a well-known hybrid functional drawback [38, 39], while the agreement between experiment and theory for SAOP and LB94 functional, specifically designed to estimate excitation energies, is somehow better. More specifically, the SAOP HOMO–LUMO ΔE (385 nm) lies very close to the absorption spectrum’s highest intensity peak (375 nm), while the LB94 one (400 nm) is at the edge of the spectrum, close to the lower energy shoulder (395 nm). As a whole, the absolute position of vertical transitions indicates that the SAOP XC functional provides the best numerical agreement between experiment and theory. Literature results pertaining to β-diketones ligands formed by a thienyl group and different polycyclic aromatic hydrocarbons (naphthyl, phenanthryl and pyrenyl) are consistent with the results herein reported [26]. Despite such a semi-quantitative agreement, we cannot be silent about the failure of the tested XC functionals in reproducing the presence of the different shoulders characterizing the absorption spectrum of L. Indeed, an accurate simulation of the excited state energies is a key factor to reproduce luminescence properties and hopes to drive the design and the synthesis of new materials.

Besides the electronic transitions, it is well known that the absorption spectrum shape is determined by different factors such as the presence of (i) diverse L geometrical isomers (such as rotamers); (ii) different L structural isomers (keto-enol tautomerism in β-diketone ligands) and (iii) vibronic progression based on the Franck–Condon principle. The shoulders’ origin in the L absorption spectrum has been then rationalized by individually considering, via DFT methods, all these factors. The absorption spectra of rotamers A, B and C have been modelled with both SAOP and LB94 functionals and, for all the rotamers, only one high-intensity TD-DFT transition (∼385 nm) is present in the UV–Vis region (see Fig. S8). The shoulders characterizing the L absorption spectrum cannot be then attributed to the presence of different L rotamers in solution and different factors need to be explored.

Even though the L NMR data are consistent with the coexistence of both enol and ketone isomers, the former species is strongly prevalent [35], perfectly in agreement with the optimized B3LYP geometry of L in its enolic form, more stable than that corresponding to the keto isomer of 6.84 kcal/mol, probably as a consequence of the higher electronic delocalisation taking place in the enol species (see in Fig. S9 the strongly localized nature of the HOMO in the ketone isomer). It has been already mentioned that a single electronic HOMO–LUMO transition at 385 nm is expected for the enol isomer on the basis of TD-DFT/SAOP calculations (see Fig. 2), while, in agreement with the literature data of Suwa et al. [40] which assigned the far UV region (∼260 nm) of the absorption spectrum of L in acetonitrile to the ketone isomer, the corresponding theoretical spectrum is characterized by several excitations around 270 nm. Thus, the ketonic form of L cannot be invoked to explain the shape of its spectrum in the region of interest (320–420 nm), which would seem to be then determined by the vibronic progression. As such, it is useful to remind that TD-DFT only considers purely electronic transitions, i.e., no vibrational components are included. The modelling of a vibrationally resolved spectrum needs the evaluation of Franck-Condon factors. Specifically, the theoretical vibrational progression has to be offset by taking into account the energy of the corresponding electronic transition, and the vibrationally resolved spectrum is obtained. Now, as the hybrid B3LYP and PBE XC functionals are not well suited to estimate excitation energies, only the LB94 and SAOP electronic transitions will be considered in the forthcoming discussion.

SAOP and LB94 vibrationally resolved absorption spectra of the rotamer A are superimposable, as they are both obtained by means of a Franck-Condon analysis carried out by adopting either the PBE or the B3LYP XC functional (see Fig. 3), the only difference being the opposite wavelength shift they need to match the recorded spectrum. Nevertheless, it has to be noted that the PBE Franck-Condon analysis is negligibly affected by the adopted number of vibrational quanta, while the B3LYP vibrational progression dramatically changes if transitions may also start from the first excited vibrational state (see Fig. 3). In even more detail, the agreement between experiment and theory is quantitative if ground state vibrationally excited levels are taken into account: the shoulders on the lower and higher wavelength side are estimated at 19 nm and 18 nm from the main peak, respectively. These two values in the experimental spectrum are 20 nm and 19 nm, respectively. As a whole, the adoption of the B3LYP functional and the inclusion of the vibronic progression not only allow a satisfactory mimicking of the absorption spectrum but also of the emission luminescence one (Fig. S11). Further details are reported in the Supplementary Information.

Fig. 3:

Normalized vibrationally resolved absorption spectra for the rotamer A of L. The calculated spectra have been redshifted by 10 nm from the HOMO → LUMO transition calculated at the TD-DFT/SAOP level. Geometries and frequencies were calculated with PBE (a) and B3LYP (b). Convoluted profiles have been obtained by using a Lorentzian broadening of 0.05 eV.

The role played by L in determining the EuL ₃ E ₂ and EuL ₃ T ₂ electronic properties and corresponding absorption spectra may be better appreciated by considering how spectral features vary upon moving from the free L to the Eu³⁺ complexes. The close similarity of the spectra (see Fig. 4a) suggests that the absorption properties are mainly determined by L, the main difference being a weak shoulder at ∼420 nm in the EuL ₃ T ₂ spectrum (see Fig. 4a).

Fig. 4:

Normalized experimental (a) and simulated (b, c) UV–Vis spectra for EuL ₃ E ₂ (b) and EuL ₃ T ₂ (c). Simulated spectra are obtained at the TD-DFT/SAOP and TD-DFT/LB94 levels by using the optimised PBE geometries. Electronic transitions of the free L have been also included for comparison..

Aimed to look into the origin of this feature, the EuL ₃ E ₂ and EuL ₃ T ₂ spectra have been then modelled. To minimize the corresponding computational cost, only electronic transitions have been considered while vibrational contributions to the spectra have been neglected. Once again, TD-DFT numerical experiments have been carried out by exploiting both LB94 and SAOP XC functionals (see Fig. 4b and c). Interestingly, at variance to L, whose 320–500 nm spectral region includes only the HOMO → LUMO transition, the spectra of both EuL ₃ E ₂ and EuL ₃ T ₂ , independently of the adopted XC functionals, are characterized by several high- and low-intensity excitations whose initial and final states are mainly generated by L-based HOMOs and LUMOs, respectively. As such, it is worthwhile to remind that the three ligands in EuL ₃ E ₂ and EuL ₃ T ₂ are not equivalent (Fig. S12). In addition, the likeness of the most intense electronic transition in EuL ₃ E ₂ and EuL ₃ T ₂ further supports the assumption that the complex spectra are scarcely affected by the Eu³⁺ presence. As far as low-intensity transitions are concerned, they both involve ligand-to-ligand charge transfer (CT) and, especially for EuL ₃ T ₂ ligand-to-metal (LMCT) transition. The LMCT minor contribution to the absorption spectra confirms that the L electronic properties are negligibly affected by the presence of the metal centre [26].

TD-DFT calculations may be also useful to look into the Eu³⁺ local environment if the attention is focused on 4f–4f transitions. Even though such analysis might be useful to thoroughly characterise the complexes, its detailed description is beyond the scope of the protocol and it is reported in the Supplementary Information.

2.3 EuL₃E₂ and EuL₃T₂ luminescence: excited states properties and intermolecular energy transfer (IET) process

The two complexes have different lifetimes and quantum yields: 60.7 μs and 0.4 % for EuL ₃ E ₂ , 170.8 μs and 5.0 %, for EuL ₃ T ₂ respectively [19]. Quite confidently, the EuL ₃ T ₂ longer lifetime may be attributed to the C₂H₅OH absence [19]. Indeed, ethanol increases the non-radiative decay rate (A_nrad) and deactivates the ⁵D₀ excited state of Eu³⁺ [19, 41, 42]. Rather than looking for the origin of the different spectroscopic properties of EuL ₃ E ₂ and EuL ₃ T ₂ , our goal is the quantitative reproduction of the experimental quantum yield PLQY and to obtain quantum efficiency, A_rad, A_nrad, the energy transfer rate (W_ET) and the back-energy transfer rate (W_BT) values by combining experimental data (lifetime and emission spectrum) and TD-DFT calculations (energies of the L lowest singlet (S₁) and triplet (T₁) excited states). In addition, the literature Eu³⁺ selected (⁵D₄, ⁵D₁ and ⁵D₀) excited states are used, whose choice has been determined by the evidence that the ⁵D₀ state may be populated by three paths: (i) the energy transfer (ET) from S₁ to ⁵D₄ and then, for non-radiative decay, to ⁵D₀; (ii) the ET from T₁ to ⁵D₁ and then, for non-radiative decay, to ⁵D₀ or (iii) a direct ET from T₁ to ⁵D₀. The weights of the different paths are related to W_ET and W_BT.

To rationalise and then quantify the different roles played by L and Eu³⁺, the energies of the L S₁ and T₁ excited states need to be known accurately. As such, the role of the triplet may be experimentally determined by substituting the Eu³⁺ ion with the Gd³⁺ one both in EuL ₃ E ₂ and EuL ₃ T ₂ : GdL ₃ E ₂ and GdL ₃ T ₂ emissions from S₁ (fluorescence) and T₁ (phosphorescence) excited states lie at ca. 23,800 cm⁻¹ and ca. 18,900 cm⁻¹, respectively [19]. As far as the theoretical modelling is concerned, a highly accurate T₁ energy is crucial because the ΔE between T₁ and the Eu³⁺ excited states ΔE(T₁ − ⁵D₀) is tightly bound to the quantum efficiency. Indeed, the T₁ state must be higher in energy than the Eu³⁺ excited state (⁵D₀) to prevent the back-ET from the Eu³⁺ emitting state to T₁, which would favour the luminescence quenching. As recently shown in a study with L = β-diketonate [26], the T₁ energy may be estimated through two distinct approaches: (i) the vertical transition (VT) approach [43] either as the lowest energy S₀ → T₁ TD-DFT vertical transition in the GS optimized geometry or as the lowest energy T₁ → S₀ TD-DFT vertical transition in the T₁ optimized geometry; (ii) the adiabatic transition – AT approach [44], as the ΔE(T₁ − S₀) with both states in their optimised geometries (see Fig. S13). AT and VT T₁ energies have been estimated for L, EuL ₃ E ₂ and EuL ₃ T ₂ as well as for LaL ₃ E ₂ and LaL ₃ T ₂ (see Table S6). According to the literature [26], the best-suited approach is the AT one and corresponding L T₁ energies will be herein used to investigate the EuL ₃ E ₂ and EuL ₃ T ₂ ET processes. A detailed analysis of VT and AT values is reported in the Supplementary Information.

To model the L − Eu³⁺ ET, both the ET and the decay rates were evaluated by employing the LUMPAC code [22], where, in addition to TD-DFT results, the experimental spectrum for the Judd-Ofelt intensity parameters (omega) and the lifetime are necessary to gain quantum yield PLQY, quantum efficiency, A_rad, A_nrad, W_ET and W_BT values. As far as Eu³⁺ is concerned, the ET has been modelled by considering only the three states ⁵D₄, ⁵D₁ and ⁵D₀ (see above), whose energy values are fixed in the code; as far as L is concerned, the S₁ and T₁ energies must be computed. Despite LUMPAC would allow their estimate, the evaluation of T₁ energies at the INDO/S-CI levels is significantly and systematically underestimated when compared with the experiment (see Table S7). Diverse geometries (PBE optimized geometry, X-ray data and Sparkle/RM1 optimized geometry) have been tested to evaluate the LUMPAC T₁ energy; its values resulted in systematically underestimated by ∼3–4000 cm⁻¹ and then a negligible energy transfer (Table S8). On the contrary, TD-DFT calculations combined with the AT approach provide T₁ values, which better agree with the experiment. On this basis, the L S₁ and T₁ energies have been estimated at the DFT/SAOP level and used as input in LUMPAC to evaluate all the parameters useful to quantitatively model the EuL ₃ E ₂ and EuL ₃ T ₂ ET mechanism. In particular, A_rad, A_nrad, quantum efficiency, and PLQY for both complexes are shown in Table 1.

The emission quantum yield PLQY value is a relevant parameter because a direct comparison with the experimental data is possible. In LUMPAC, PLQY value is the ratio between the emitted and absorbed light intensities (see Eq. (3)):

where η D 0 5 and η S 0 are the ⁵D₀ and S₀ level populations, while φ correspond to absorption rate. The ⁵D₀ level population depends on the non-radiative decay rate A_nrad, which can be obtained from A_rad and the experimental lifetime (τ) as A_rad + A_nrad = τ⁻¹. The radiative decay rate (A_rad) directly depends on the Judd-Ofelt intensity parameters. Thus, the theoretical PLQY is related to the experimental lifetime ad emission spectrum.

The inspection of data reported in Table 1 reveals the substantial difference between A_rad and A_nrad values for both complexes. Moreover, the EuL ₃ E ₂ A_nrad value (15814 s⁻¹) is larger than the A_rad one (661 s⁻¹), while the opposite is true for EuL ₃ T ₂ (4712 s⁻¹ and 1143 s⁻¹ for A_nrad and A_rad, respectively). Such a difference can be explained by referring, in the former complex, to the presence of the O–H oscillators, which accounts for the lower PLQY and quantum efficiency. Theoretical outcomes herein reported allow us to gain further details on the radiative and non-radiative process: the EuL ₃ E ₂ A_rad value is ∼½ the EuL ₃ T ₂ one, while the EuL ₃ E ₂ A_nrad is ∼3 the EuL ₃ T ₂ one. This result suggests that ancillary ligands such as EtOH or TPPO influence differently radiative and non-radiative processes: they slightly affect the former, while they may significantly affect the latter. Further information about A _rad may be gained by referring to the contribution percentage of each Eu^{3+ 5}D₀ → ⁷F_J transition; in fact, not all transitions equally contribute to the radiative process. In this regard, it can be useful to mention that A_rad and A_nrad have been calculated in the past for similar Eu³⁺ complexes containing TPPO by using LUMPAC and literature results show similar trends [45, 46]. In the EuL ₃ E ₂ (EuL ₃ T ₂ ) percentages to the radiative decay rate are ⁵D₀ → ⁷F₁ = 7.45 (4.31), ⁵D₀ → ⁷F₂ = 73.13 (88.52), ⁵D₀ → ⁷F₃ = 0.00 (0.00), ⁵D₀ → ⁷F₄ = 19.39 (7.14), ⁵D₀ → ⁷F₅ = 0.00 (0.00), and ⁵D₀ → ⁷F₆ = 0.03 (0.03). These values show that, in both complexes, the ⁵D₀ → ⁷F₂ transition has the largest percentage thus providing the major contribution to the radiative process, even if with sizable different percentages. Finally, the populations of the S₀, ⁵D₀ and T₁ states have been calculated for both complexes: η S 0 , η D 0 5 , η T 1 are 0.844, 0.086 and 0.069, respectively, for EuL ₃ E ₂ and 0.721, 0.220 and 0.059, respectively, for EuL ₃ T ₂ . It is interesting to highlight that the population on S₀ and T₁ levels are quite similar for the two complexes, while the largest variation is on the ⁵D₀ level. This large variation in the population of the emitter level difference justifies the larger calculated quantum yield PLQY value for EuL ₃ T ₂ (3.48%) with respect to EuL ₃ E ₂ (0.67 %). These values are very close to the experimental ones equal to 5.0% and 0.4%, for EuL ₃ T ₂ and EuL ₃ E ₂ , respectively.

As far as W_ET and W_BT are concerned, they have been calculated within the Malta’s model [17, 47] by using LUMPAC code (see Table 2). In more detail, W_ET and W_BT are obtained for the most relevant processes involving the ligand and the Eu³⁺ sides, S₁ ↔ ⁵D₄, T₁ ↔ ⁵D₁ and T₁ ↔ ⁵D₀ and they are presented in Table 2 and in Figs. S14 and S15.

The position of the excited states involved in the ET mechanism makes W_ET and W_BT very similar in both EuL ₃ E ₂ and EuL ₃ T ₂ . The first possible path to populate the ⁵D₀ state is the ET from S₁ to ⁵D₄. The S₁ → ⁵D₄ W_ET is 10⁴ s⁻¹, definitely slower than the S₁ → T₁ inter-system crossing (rate = 10⁸ s⁻¹) [48], which cancels the former path weight. Such a peculiar behaviour is determined by the relative position of the ⁵D₄ state, higher in energy than the S₁ one by ∼2000 cm⁻¹. The T₁ state is populated at the expense of the ⁵D₄ one and the main ET paths are from L T₁ to Eu³⁺ lower excited states (⁵D₁ and ⁵D₀). A second possible alternative path to populate the ⁵D₀ state is the ET from T₁ to ⁵D₁. For both complexes, the T₁ → ⁵D₁ W_ET is 10⁴ s⁻¹, while W_BT amounts to 10⁶ s⁻¹. The latter value is of the same order of magnitude as the ⁵D₁ → ⁵D₀ non-radiative rate; thus, back-ET and the non-radiative decay are competitive processes and the ⁵D₁ state may contribute to populating the ⁵D₀ one. The third and final path corresponds to the direct T₁ → ⁵D₀ ET. For both complexes, the corresponding W_ET = 10⁴ s⁻¹, while W_BT = 10³ s⁻¹. The direct ET seems thus to be the process most contributing to the ⁵D₀ state population. Before going on, it can be of some interest to point out that the T₁ → ⁵D₀ W_ET is only twice the T₁ → ⁵D₁ one. This is due to the peculiar T₁ relative position, whose energy is almost in between the ⁵D₀ and ⁵D₁ ones (see Figs. S14 and S15). Based on the theoretical luminescence parameters, we propose the following EuL ₃ E ₂ and EuL ₃ T ₂ ET diagram (Fig. 5). Dashed and full lines are associated with nonradiative and radiative paths, respectively, while curve dashed lines account for L → Eu³⁺ ET or Eu³⁺→ L back-ET.

Fig. 5:

EuL ₃ T ₂ and EuL ₃ E ₂ Jablonski energy-level diagram with most relevant ET channels. Theoretical photophysical properties (see also Tables 1 and 2) have been obtained by exploiting the LUMPAC code. Red and black arrows account for L → Eu³⁺ ET and Eu³⁺→ L back-ET, respectively. W_ET and W_BT values are reported in Table 2.

In general, an efficient ligand → metal ET implies in Eu³⁺ complexes implies 2500 cm⁻¹ < ΔE(T₁ − ⁵D₀) <4000 cm⁻¹ [29]. Thus, in the presence of β-diketonate ligands, the ET is dominated by the exchange Dexter mechanism and W_ET is in the order of 10⁸ s⁻¹ [41, 49]. Moreover, in several diketonate complexes, E(T₁) > E(⁵D₀) and E(T₁) > E(⁵D₁) so that the T₁ → ⁵D₁ ET rate becomes the largest one [41]. Conversely, the T₁ state has in the present case a peculiar position relative to the Eu^{3+ 5}D₀ and ⁵D₁ low-lying excited states. Indeed, E(T₁) > E(⁵D₀) by only ∼1800 cm⁻¹ but E(T₁) < E(⁵D₁) by ∼1100 cm⁻¹. This evidence implies: (i) T₁ → ⁵D₀ and T₁ → ⁵D₁ W_ET are rather low (∼10⁴ s⁻¹) and (ii) the T₁ → ⁵D₁ W_ET is not the largest. Even though a further consequence of the peculiar T₁ relative energy is the low quantum yields, it has to be underlined that EuL ₃ E ₂ and EuL ₃ T ₂ are extremely sensitive to temperature variations and thus making them useful for technological applications such as molecular thermometers.