Ultrafast dynamics and excited-state trapping in [3.3]paracyclophane

Abstract

Paracyclophanes are rigid, three-dimensional frameworks in which two benzene rings are held in a parallel, stacked arrangement by short aliphatic linkers. Their derivatives display pronounced through-space π–π interactions that strongly influence their excited-state behavior. Here, we investigate the ultrafast excited-state dynamics of [3.3]paracyclophane ([3.3]PCP) using time-dependent density functional theory combined with nonadiabatic molecular dynamics via surface hopping. The 2 ps simulations provide a detailed picture of how electronic and nuclear motions evolve in concert after photoexcitation. Following excitation to the weak S₃ band shoulder, [3.3]PCP rapidly relaxes into S₁, where it becomes kinetically trapped. Analysis of the fragment-based transition density matrix reveals a concurrent transformation of the electronic structure—from an excitonic-resonance state with minor charge-transfer (CT) character (≈0.2) to a mixed excitonic/charge-resonance regime with CT ≈ 0.5. This evolution is accompanied by a structural contraction of the π-stack, as the two benzene rings approach each other, activating an inter-ring breathing motion that governs the subsequent dynamics. Because [3.3]PCP combines rapid nonadiabatic relaxation with long-lived excited-state trapping, it serves as a particularly demanding benchmark for trajectory-based dynamics methods, offering both a mechanistic picture of excimer formation and a stringent test of their current capabilities.

---

1. Introduction

Designing efficient organic electronic materials often comes down to a delicate balance: molecules must couple strongly through their π-systems, yet avoid the conformational disorder that undermines transport.¹ Rigid frameworks are one way to achieve this balance.² Paracyclophanes (PCPs) are a striking example (Fig. 1). Their aliphatic linkers force two benzene rings into close parallel alignment, creating short through-space π–π contacts while interrupting through-bond conjugation between the rings.³ The result is a scaffold with distinctive excited-state interactions arising from enforced π–π proximity.⁴ Derivatives of [3.3]PCP have already found use as donors and acceptors in thermally activated delayed fluorescence emitters,^5–7 as rigid spacers in supramolecular assemblies,⁸ as radical stabilizers,⁹ and as elements in molecular junctions where quantum interference boosts conductance.¹⁰

Fig. 1 Molecular structure of [3.3]paracyclophane ([3.3]PCP).

PCPs are also proving valuable as model systems. In parallel with the present study, Haggag and co-workers combined ultrafast spectroscopy with high-level calculations to show that [2.2]PCP displays excimer-like behavior:¹¹ excitation triggers an inter-ring breathing motion that modulates excimer absorption on the femtosecond scale. This result illustrates how the enforced geometry of cyclophanes provides a controlled setting for probing fundamental excited-state processes.

Here, we focus on [3.3]PCP and investigate its excited-state dynamics using nonadiabatic trajectory surface-hopping simulations. This approach allows us to follow how the electronic character evolves after photoexcitation on the picosecond timescale. In contrast to Haggag et al.'s work,¹¹ which used femtosecond spectroscopy to capture coherent vibrational modulation, our study resolves the nonadiabatic relaxation and long-lived interplay between excitonic and charge-resonance states. Because [3.3]PCP undergoes rapid population transfer through higher-state conical intersections but subsequently becomes trapped in a regime of sustained S₁/S₂ exchange, it provides a particularly demanding test case for current nonadiabatic dynamics methods.

To this end, we employ trajectory surface hopping at the TDDFT level, which captures the picosecond-scale relaxation following photoexcitation and traces how electronic populations respond to nuclear motion. To interpret these trajectories, we use a fragment-based transition-density framework that provides a direct view of the excitonic and charge-resonance character along the relaxation pathway.¹² This combined strategy reveals how structural rearrangements of the cyclophane scaffold mediate the transformation from an initial local excitation into long-lived states of mixed excitonic and charge-resonance character, while simultaneously providing a demanding test case for the reliability of widely used excited-state methods. Although surface hopping is not expected to give a complete description of long-lived coherence,¹³ it serves as a critical baseline: a first step toward assessing how current approaches perform in systems explicitly designed to challenge them.

2. Computational details

Ground-state calculations were performed using DFT with two functionals, CAM-B3LYP¹⁴ and ωB97X.¹⁵ Excited-state calculations were done using linear-response TDDFT with these functionals. These calculations were performed using the quantum chemistry software Turbomole V7.6.¹⁶ All calculations were conducted with the basis set def-SV(P).¹⁷ The resolution of identity (RI) approximation was also employed in all calculations, as it is implemented in Turbomole, and it can significantly reduce computational cost. Long-range dispersion interactions (van der Waals forces) were accounted for using the D4 dispersion correction model¹⁸ throughout the calculations. Cartesian coordinates of the stationary points are given in the SI.

Additional calculations were performed using TDDFT with the Tamm–Dancoff approximation (TDDFT-TDA) using the CAM-B3LYP functional, DFT-based multireference configuration interactions (DFT/MRCI),^19,20 algebraic diagrammatic construction to the second order (ADC(2)),²¹ coupled cluster with approximated second order (CC2),²¹ and rotated multistate complete active space to the second order perturbation theory (RMS-CASPT2).²² DFT/MRCI was run with the Hamiltonian and parameters proposed in ref. 23. The energy selection criterion was set to 1 au. RMS-CASPT2 vertical excitation energies were computed using a CAS(12,12) (complete active space) reference, including three π and three π* orbitals for each of the benzene rings. The CASSCF orbitals were obtained from a state-averaged CASSCF over 11 singlet states, and the perturbative correction was evaluated with the RMS-CASPT2 variant (rotated multistate CASPT2). An IPEA shift of 0.0 and an imaginary level shift of 0.2 a.u. were used to control intruder-state effects. All TDDFT, TDDFT-TDA, DFT/MRCI, ADC(2), and CC2 calculations employed the SV(P) basis set. The RMS-CASPT2 calculations employed the cc-pVDZ basis set.

TDDFT and TDDFT-TDA, ADC(2), and CC2 calculations were performed using Turbomole. DFT/MRCI was computed with the DFT/MRCI program^19,20 using DFT ground-state calculations with Turbomole. For RMS-CASPT2, we used OpenMolcas.²⁴ We searched for S₁/S₀ interstate crossings using the penalty-function approach in an in-house modified version of CIOpt.²⁵

To initiate nonadiabatic dynamics, an ensemble of 1800 initial conditions was generated using a harmonic Wigner distribution centered around the optimized ground-state (S₀) minimum geometry, computed with the B3LYP functional²⁶ and the def-SV(P) basis set. For each sampled geometry, vertical excitation energies and oscillator strengths were calculated using TDDFT/CAM-B3LYP+D4/def-SV(P). This information was used to simulate the absorption spectrum using the nuclear ensemble approach.²⁷ A Gaussian broadening of 0.05 eV was applied. The initial conditions and spectrum were calculated using the CS module of the Newton-X platform,²⁸ interfaced with Turbomole.

Nonadiabatic dynamics were performed using the Fewest Switches Surface Hopping (FSSH)²⁹ with the Granucci-Persico simplified decay of mixing decoherence corrections (DC; α = 0.1 au).³⁰ The dynamics were run for up to 2 ps, using a classical time step of 0.5 fs, resulting in up to 4000 single-point calculations per trajectory. Quantum equations were integrated with a time step of 0.025 fs, using interpolated electronic properties. The time-derivative nonadiabatic couplings required for FSSH were computed using the time-dependent Baeck-An (TD-BA) approximation,³¹ which extends the curvature-based model of Baeck and An³² to a time-dependent form. Different algorithms for velocity rescaling after a hop were employed: velocity rescaling in the gradient-difference vector direction (g approach), in the momentum direction (p approach), and in the momentum direction with reduced kinetic energy reservoir (p-redKE approach).³³ A total of 100 trajectories were simulated for each approach. The g approach, implemented in Newton-X for this work, incurs almost no additional cost, as it only requires evaluating additional gradients when hopping events occur, rather than at every step. In frustrated hopping, momentum was conserved in the original direction. DC-FSSH was performed using the NS module of the Newton-X platform,²⁸ interfaced with Turbomole.

The ring conformations were analyzed with the Cremer-Pople-Boeyens approach. Cremer-Pople analysis is a method for quantifying and describing the puckering of ring molecules.³⁴ It projects the ring's conformation onto a set of parameters, which tells how the ring is distorted. The puckering of six-membered rings is described in terms of three parameters: the puckering amplitude Q, and the angular parameters θ and ϕ, which can be mapped on Boeyens’ qualitative ring distortion types:^35,36 boat (^ijB or B_ij), twisted boat (ⁱT_j), screw-boat (ⁱS_j), chair (ⁱC or C_i), half-chair (ⁱH_j), and envelope (ⁱE or E_i) involving atoms i and j. Subscripts indicate a puckering of atoms below the plane, while superscripts denote puckering above the plane. Structures with Q smaller than 0.05 Å were considered planar.

To monitor the evolution of the excited-state electronic character along the trajectories, we analyzed the one-electron transition density matrix (1TDM) between the ground state |Ψ₀〉 and the current excited state |Ψ_α〉. The 1TDM is defined as¹²

where r and s denote orbital indices and â^†, â are the usual creation/annihilation operators.

Following Plasser and Lischka,¹² we constructed the charge-transfer number matrix Ω^(0α), whose elements Ω^(0α)_AB quantify the weight of configurations where the hole is located on fragment (A) and the electron on fragment (B). In an orthonormal localized-orbital (LO) basis, this quantity can be written as

In practice, Ω^(0α) was computed with TheoDORE using the AO-based generalization that accounts for AO overlap (Mulliken/Mayer-type partitioning).³⁷ The total norm is

which is expected to be close to unity for states dominated by one-electron transitions.

From Ω, we evaluated two scalar descriptors: the charge-transfer character (CT) and the participation ratio (PR). The CT value is the total off-diagonal weight,

so that CT ≈ 0 corresponds to predominantly local/Frenkel-type character (electron and hole on the same fragment), while CT ≈ 1 corresponds to predominantly charge-separated character.¹² PR quantifies delocalization. It is defined as

PR ≈ 1 indicates localization on a single fragment, whereas PR ≈ 2 indicates delocalization over two fragments; values >2 indicate participation of additional fragments.

In the present work, CT and PR were computed by partitioning [3.3]PCP into four fragments as discussed in Section 3.4. The state classification according to these descriptors is also discussed in that section.

All statistical analyses, including the Cremer–Pople–Boeyens approach, were done with ULaMDyn.³⁸

3. Results and discussion

3.1. Static results and benchmark calculations

We began by optimizing [3.3]PCP at its ground state. Normal modes were calculated and checked to ensure that there were no imaginary frequencies, indicating a minimum. The benzene rings are nearly planar with the inter-ring distance of 3.25 Å. The IR spectrum is provided in the supplementary information (SI) for reference. Vertical excitation energies for the first 10 excited states for [3.3]PCP were calculated using TDDFT (CAM-B3LYP), TDDFT (ωB97X), TDDFT-TDA (CAM-B3LYP), DFT/MRCI (on the CAM-B3LYP geometry), ADC(2), CC2, and RMS-CASPT2 (on the CAM-B3LYP geometry) and, as shown in Table 1.

Table 1 Vertical excitation energies (ΔE in eV) and oscillator strengths (in parentheses) for the first 10 singlet excited states for [3.3]PCP in the S₀ minimum computed with different electronic structure methods

State	TDDFT (CAM-B3LYP)	TDDFT-TDA (CAM-B3LYP)	TDDFT (ωB97X)	DFT/MRCI	ADC(2)	CC2	RMS-CASPT2 (12,12)
	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)	ΔE (eV) (Osc. str.)
S1	4.68 (0.000)	4.73 (0.000)	4.82 (0.000)	4.32 (0.000)	4.52 (0.000)	4.50 (0.000)	4.14 (0.000)
S2	5.06 (0.000)	5.21 (0.000)	5.28 (0.000)	4.90 (0.008)	5.10 (0.002)	5.07 (0.002)	4.64 (0.001)
S3	5.41 (0.006)	5.45 (0.005)	5.44 (0.006)	4.96 (0.000)	5.32 (0.000)	5.31 (0.000)	4.65 (0.000)
S4	5.61 (0.000)	5.83 (0.000)	5.88 (0.000)	5.42 (0.000)	5.79 (0.000)	5.78 (0.000)	5.00 (0.000)
S5	5.77 (0.000)	5.97 (0.000)	6.02 (0.000)	5.56 (0.000)	5.99 (0.000)	5.97 (0.000)	5.16 (0.000)
S6	6.12 (0.122)	6.31 (0.123)	6.17 (0.110)	5.92 (0.160)	6.23 (0.203)	6.20 (0.187)	7.02 (0.000)
S7	6.42 (0.001)	6.41 (0.003)	7.05 (0.015)	6.20 (0.001)	6.67 (0.016)	6.60 (0.018)	7.20 (0.000)
S8	6.67 (0.001)	6.63 (0.002)	7.25 (0.103)	6.43 (0.004)	6.81 (0.059)	6.75 (0.055)	7.26 (0.000)
S9	6.93 (0.013)	6.89 (0.011)	7.29 (0.551)	6.63 (0.041)	7.04 (0.132)	6.99 (0.124)	7.33 (0.000)
S10	7.10 (0.061)	7.09 (0.008)	7.38 (1.340)	6.89 (0.000)	7.17 (0.052)	7.09 (0.000)	7.51 (0.000)

---

Across all methods, the lowest singlet excitations remain essentially dark at the S₀ minimum (Table 1): up to S₅, oscillator strengths are negligible (typically f ≲ 10⁻²), whereas a much more intense band emerges from S₆ (f ≈ 0.1–0.2 in TDDFT, DFT/MRCI, ADC(2), and CC2). This trend is consistent with an approximate C_2h description of the Franck–Condon geometry, for which the lowest states are predominantly gerade-like and therefore dipole-forbidden. At the same time, weak intensity appears only once a B_u-like state (dipole-allowed along z, the C₂ axis) is reached. In this respect, TD-CAM-B3LYP, TD-CAM-B3LYP-TDA, TD-ωB97X, and DFT/MRCI place the first weakly allowed state at S₃, whereas ADC(2), CC2, and RMS-CASPT2 place it at S₂. Notably, the corresponding oscillator strengths remain small in all cases, and the overall energetic ordering of the low-lying manifold is similar.

Excitation energies are systematically lower with DFT/MRCI and higher with TD-ωB97X relative to TD-CAM-B3LYP, with ADC(2) and CC2 lying in between. TDDFT-TDA results closely follow full TDDFT values for both energies and oscillator strengths, except for minor shifts (≤0.2 eV). RMS-CASPT2 reproduces the existence of a dark low-energy region; however, within the present protocol, it predicts a noticeably different energy scale for higher excitations (from S₆ upward), likely reflecting the combined sensitivity to the reference geometry and the chosen multireference/PT2 setup (active space and state averaging). The significant spread in f values for higher states (e.g., S₉–S₁₀) suggests varying descriptions of transition character among methods, with TD-ωB97X predicting particularly intense transitions for S₉ and S₁₀. Analysis of the DFT/MRCI wave functions shows that all states S₁–S₁₀ are dominated by single excitations, with single-excitation weights between 89 and 94%.

Taken together, the cross-method consistency in the low-energy manifold (S₁–S₅) at the Franck–Condon region supports using TD-CAM-B3LYP as a practical working level for the nonadiabatic dynamics—primarily because it is computationally feasible for this system—while acknowledging that the resulting timescales should be interpreted as method-dependent rather than strictly benchmarked values.

In the S₁ minimum, the benzene rings remain nearly planar, but the inter-ring distance reduces to 2.94 Å. In this case, in the approximate C_2h description, the S₁ state has B_g symmetry. We have also located an S₁/S₀ intersection geometry, which is discussed later (Section 3.3). Energy information about the minima and the intersection point is given in Table 2. Vertical absorption occurs at 265 nm while vertical emission occurs at 322 nm, corresponding to a 57 nm Stokes shift.

Table 2 Energies of [3.3]PCP at different geometries computed with TDDFT (CAM-B3LYP). ΔE_v is the vertical energy gap. Oscillator strengths are given in parentheses. For the S₀ minimum, the results are the same as in Table 1

Geometry	E(S0) (eV)	ΔEv(S1–S0) (eV)	ΔEv(S2–S0) (eV)
S0 minimum	0.00	4.68 (0.000)	5.06 (0.000)
S1 minimum	0.46	3.85 (0.000)	4.19 (0.000)
S1/S0 intersection	7.47	0.03	2.44

---

3.2. Surface hopping dynamics

The absorption spectrum of [3.3]PCP was simulated using the nuclear ensemble approach, considering excitations into the first nine singlet excited-state transitions (S₁–S₉). The spectrum displays two main bands (Fig. 2). The dominant band, centered at ∼6.1 eV, originates from the bright S₆ transition, as evident from the spectrum that considers excitations into this state exclusively. A weaker band appears at ∼6.6 eV and is attributed to higher excited states (S₇–S₉). In addition, the main 6.1 eV band exhibits a shoulder on its low-energy side at ∼5.5 eV, arising from the S₃ state.

Fig. 2 Absorption cross section (Ų per molecule) of [3.3]PCP considering the first nine excited states. The peaks corresponding to S₃ and S₆ only are also shown.

The initial conditions were selected from the spectral shoulder between 5.00 and 5.50 eV, centered at 5.25 eV (±0.25 eV). To ensure physical relevance, only initial conditions corresponding to excited states with significant oscillator strengths were retained.²⁸ The dynamics simulations were initiated in S₃. Each trajectory was allowed to evolve nonadiabatically between the first five singlet excited states (S₁–S₅). Explicit hops to the ground state were not evaluated. However, when the S₁–S₀ energy gap became smaller than 0.2 eV, transition to the ground state was assumed to occur. With this information, it is possible to reconstruct the S₀ occupation.

After excitation, the S₃ population quickly transfers to S₂ within ∼50 fs (Fig. 3). S₂ immediately transfers most of the population to S₁, which remains trapped in this state until the end of the dynamics, with approximately 10% to 80% of the population occupying states S₂ and S₁, respectively. There is an indication of a slow transfer to the ground state on a much longer time scale. We return to this point later. The characteristic time constants for diverse processes are compiled in Table 3.

Fig. 3 Mean state occupation over time from surface hopping dynamics using momentum rescaling along the gradient difference vector (g). On the left side, a zoomed-in plot of the mean occupation for the first 200 fs, and on the right side, the whole plot for 2 ps.

Table 3 Time constants for different processes during [3.3]PCP dynamics. The S₀ population time constant can be interpreted as an effective global excited-state lifetime within this surface-hopping model

Feature	Time constant (ps)
S3 depopulation	∼0.05
S0 population	19 [12–23]
Mixed ER + CR formation	∼0.04

---

Proper velocity adjustment after electronic transitions is critical to avoid artifacts in surface hopping, particularly size-extensivity issues.³⁹ The results in Fig. 3 were computed rescaling velocities in the gradient difference direction (g), which is recommended when nonadiabatic coupling vectors are unavailable.³³ We also run dynamics with velocity rescaling in the momentum direction (p) and in the momentum direction with reduced kinetic energy (p-redKE). These results are presented in the SI S2.

The g and p-redKE schemes yield nearly identical population dynamics: rapid relaxation into S₁ followed by long-lived trapping, with only a modest, slow rise in the ground-state occupation. In contrast, the simple momentum (p) rescaling yields a qualitatively different outcome, with an artificial excess of backhopping, which distorts the entire population dynamics. This behavior arises from the vast kinetic-energy reservoir made available in the p scheme, consistent with our previous analysis.³³ Throughout the rest of the paper, we exclusively focus on the g-rescaled dynamics.

Before discussing these results, we want to comment on the applicability of the TD-BA approximation to [3.3]PCP. The original BA model was derived for a one-dimensional, locally two-state problem,³² and the S₁–S₅ manifold in [3.3]PCP is indeed relatively dense, spanning only about 1.1 eV at the Franck–Condon geometry (Table 1). This, however, does not imply the presence of genuine three-state (or higher) conical intersections along the relaxation pathway, which are nongeneric and typically require additional symmetry or fine-tuning.⁴⁰ In DC-FSSH, nonadiabatic events are treated pairwise between the active state and a single target state, consistent with the TD-BA formalism, which reconstructs the effective coupling for each pair from the local time dependence of their energy gap. Previous works showed that, in such a pairwise, time-local regime, TD-BA yields a qualitatively correct description and semi-quantitative agreement with the exact nonadiabatic couplings for multistate dynamics.^31,41 Thus, we regard TD-BA as introducing an additional, but moderate, source of uncertainty in the absolute values of the time constants reported here, while leaving the qualitative mechanism unchanged, as recently discussed in ref. 42.

To model the kinetics of ground-state repopulation, the time evolution of S₀ occupation was fit to a single-exponential growth function with a delayed onset. The chosen model accounts for the fact that the S₀ occupation remains flat for an initial period, followed by a smooth, asymptotic increase. The fitting function is defined as:

where P(t) is the S₀ occupation, and t₀ is the onset time at which the occupation begins to rise, and τ is the time constant characterizing the rate of occupation growth. This same time constant can be interpreted as a global excited state lifetime. This model assumes complete relaxation to the ground state (final occupation = 1.0) and mimics the behavior observed in the surface-hopping trajectories: S₀ remains essentially unpopulated during an initial induction period, and once trajectories reach the S₁/S₀ crossing region, the S₀ population grows approximately exponentially. In this context, the time constant τ represents an effective lifetime of the excited-state ensemble.

Fitting the S₀ occupation with eqn (3) yielded an onset time of t₀ = 0.13 ps and a lifetime of τ = 19 ps (Fig. 4), indicating that ground-state recovery is a relatively slow process. Because the simulations were limited to 2 ps, this value represents an extrapolation rather than a directly sampled decay. To estimate its reliability, we applied bootstrap analysis to compute the 95% confidence interval, yielding bounds of 12–23 ps. This timescale should therefore be viewed as an order-of-magnitude estimate of the S₁ → S₀ decay for an isolated [3.3]PCP molecule within the present TDDFT/DC-FSSH framework and the adopted small-gap hopping criterion, rather than as a high-precision prediction of the experimental lifetime.

Fig. 4 Exponential fit of the mean ground-state occupation S₀ extracted from surface-hopping dynamics. The dashed orange curve represents the fitted model with onset time t₀ = 0.13 ps and lifetime τ = 19 ps. The shaded area shows the 95% confidence interval obtained from the bootstrap analysis, corresponding to a lifetime uncertainty of 12–23 ps. Because the dynamics were propagated for only 2 ps, the fitted time constant represents an extrapolation beyond the simulated time window.

The long time constant for S₀ population growth (or excited-state depopulation) reflects the kinetically hindered nature of S₁ → S₀ internal conversion, indicating that while nonadiabatic relaxation from initially excited S₃ state occurs within 0.1 ps, ground-state recovery is the slowest step in the relaxation cascade. This is consistent with a bottlenecked relaxation pathway dominated by S₁ trapping, requiring tens of picoseconds to fully repopulate the ground state.

3.3. Structural analysis

To investigate how the [3.3]PCP molecule behaves throughout its dynamics, we analyzed its ring conformation using the Cremer–Pople–Boeyens approach, as shown in Fig. 5.

Fig. 5 (left side) Division of [3.3]PCP for Cremer–Pople–Boeyens analysis. (right side) Heatmap showing the data of puckering amplitude (Q) for both (upper and lower) rings.

The Cremer–Pople–Boeyens analysis was conducted independently for each ring at each time step throughout all trajectories. The heatmap of puckering amplitudes (Q) shown in Fig. 5 reveals that the mean amplitude for each ring is approximately 0.18 Å, with values typically ranging between 0.05 and 0.25 Å. The density is concentrated along the Q_upper ≈ Q_lower diagonal, and structures where one ring is strongly puckered while the other is nearly planar are rare. This pattern indicates a positive correlation between the puckering amplitudes of the two rings.

For clarity, we adopt an ortho/meta/para notation for puckering atoms relative to the para-substitution pattern of benzene, grouping all symmetry equivalent conformations into the same class; for instance, B_2,5 is denoted B_p,p, ³S₂ as ^mS_p, and ⁶T₂ as oT_p.

Fig. 6 highlights the most significant conformations observed during the excited-state dynamics. B_p,p, with C2 and C5 puckered toward the opposite ring, dominates. Other significant conformers include the twist-boat (^mT_p) and envelope-like forms (E_p), both of which have atoms C2 or C5 puckered toward the opposite ring. The presence of multiple distinct puckering conformations involving C2 or C5 puckering suggests a restricted conformational flexibility during the dynamics. In fact, the dominant dimer conformation is ^p,pB–^mT_p (31k occurrences out of 387k structures analyzed), followed by ^p,pB–^mS_p (24k), ^p,pB–B_p,p (23k), and ^p,pB–E_p (14k). Fully planar conformers have a meager 529 occurrences.

Fig. 6 Histogram of the Cremer–Pople–Boeyens analysis for both rings for all trajectories over all times. In these conformations, the para (p) carbon atom is always puckered toward the opposite ring.

During the excited-state dynamics, the mean distance between the phenyl rings contracts from 3.25 Å, the equilibrium value in S₀, to 3.04 Å, close to the S₁ minimum at 2.94 Å, within about 200 fs. It then oscillates with a frequency of 172 cm⁻¹ (period ≈ 194 fs; Fig. 7). For roughly a picosecond, the individual trajectories oscillate in phase; subsequently, small differences in their frequencies cause them to drift out of phase, so that the ensemble-averaged inter-ring distance loses its apparent oscillations.

Fig. 7 Mean inter-ring distance as a function of time. The shaded region marks the standard deviation.

Normal-mode analysis at the S₁ minimum shows that this vibration is the inter-ring breathing mode, calculated at 202 cm⁻¹. In the ground state, the same mode lies at 169 cm⁻¹, indicating that excitation tightens the inter-ring potential—precisely what one expects when an excimer bond forms between the π systems. The lower frequency observed in the dynamics (about 0.85 of the harmonic value) reflects anharmonic softening and coupling to low-frequency torsional motions during the large-amplitude excursion along this coordinate.

This breathing motion is the same one identified by Haggag et al. in [2.2]PCP, where coherent oscillations near 233 cm⁻¹ on the S₁ surface modulate the excimer absorption.^11,43 The smaller value found here for [3.3]PCP (172 cm⁻¹) naturally fits with its longer, more flexible bridges and weaker inter-ring coupling, with the excimer spring relaxed by an additional methylene unit.

The slow radiationless S₀ recovery is mediated by an S₁/S₀ state intersection characterized by a ^p,pB–B_p,p configuration, corresponding to a boat distortion where both rings pucker toward each other. Its geometry is illustrated in SI S3, and Cartesian coordinates are in S4. This intersection, lying 3.2 eV above the S₁ minimum (Table 2), is reached in an uphill topography, which is unfavorable for internal conversion. Note that this minimum-energy crossing point was optimized under the strict condition ΔE(S₁–S₀) ≤ 0.03 eV and therefore marks the bottom of the crossing seam. Along the inter-ring breathing coordinate, configurations with small but finite S₁–S₀ gaps (ΔE ≲ 0.2 eV) lie significantly lower in energy than this more strict crossing point and can be transiently reached with the available vibrational excess. The rare S₁ → S₀ hops observed in the dynamics are associated with visits to these small-gap regions rather than to the strict crossing itself. SI S3 shows an example of a trajectory reaching the ΔE ≲ 0.2 eV crossing region.

The linear interpolation path between the ground-state minimum and the intersection geometry is provided in the SI S3. Note that our dynamics is for an isolated molecule. As such, it does not account for the dissipation of vibrational energy to the environment. If this effect were accounted for, the radiationless ground state recovery would likely be significantly reduced or even disappear entirely.

3.4. Fragment-based transition density matrix analysis

To gain further insight into the [3.3]PCP excited-state dynamics, we investigated how electronic density and excitonic character evolve by analyzing the 1-electron transition density (1TDM). We analyzed the 1TDM distribution into four fragments: fragments 1 and 2 are the two benzene rings, and fragments 3 and 4 are the side chains, as shown in Fig. 8.

Fig. 8 Split of [3.3]PCP for the fragment-based transition density matrix analysis.

We chose 28 points (10 with equal intervals in the first 200 fs, as most state transitions occur within this time frame, and 18 with equal intervals from the next 1800 fs) during the dynamics to investigate how the charge character of the excited states changes. Two 1TDM descriptors—CT (eqn (1)) and PR (eqn (2))—were monitored at these time steps for all trajectories. They were used to characterize the excited-state electronic configuration. Fig. 9 schematically depicts five different electronic states that can be expected to be formed in a similar stacked π system.^37,44

Fig. 9 Schematic representation of the different pure and mixed excited states that can be formed in [3.3]PCP with their expected charge transfer Ω-matrix and the respective values of the charge transfer (CT) and participation ratio (PR) descriptors.

Pure local and excitonic states are characterized by the absence of charge transfer, i.e., the off-diagonal elements of the Ω-matrix are close to zero. While local excitations (LE) are localized on a single fragment, which for [3.3]PCP is the benzene ring, excitonic resonance (ER) occurs when the excitation is delocalized over two or more fragments, which is distinguished through their PR. Analogously, two types of pure charge transfer states can be formed, characterized by a full charge transfer transition (CT = 1), where the diagonal elements of the Ω-matrix are zero. A charge transfer state (CTS) occurs when the excitation occurs from one fragment to another, so the hole and electron are each localized on a single, different fragment, resulting in PR = 1. Charge resonance (CR) is achieved when charge transfer from each fragment to the other occurs, resulting in a net zero charge transfer. Finally, a mixed regime between excitonic and charge resonance (Mixed ER + CR) can be formed. This will be characterized through a CT value of 0.5.

As illustrated in Fig. 9, all five types of excited states can be clearly differentiated by their Ω-matrix and the PR and CT descriptors. We defined cut-off thresholds for the descriptors to create a classification scheme for the excited states. For PR, a threshold of 1.5 was used; states with PR < 1.5 were classified as local states (either LE or CT states), while ER, CR, and mixed ER + CR were defined for PR ≥ 1.5. Localized states, LE and ER, were defined for CT < 0.3. Pure CT and CR states were defined for CT ≥ 0.7, while mixed ER + CR were defined in between, for 0.3 ≤ CT < 0.7.

As illustrated in Fig. 10, the majority of states along the trajectories can be classified as mixed ER + CR. A significant clustering of the states occurs centred around CT = 0.5 and between PR of 2.0 and 2.75. The larger values of PR, many exceeding 2.0, indicate that the excitation delocalises not only over the two benzene rings but also involves the side chains, which are treated as separate fragments. There are fewer points where the molecule is in an ER state, and even fewer in an LE state. No CR or CT states are observed at any point along the dynamics.

Fig. 10 Distribution of points during the dynamics in terms of the PR and CT descriptors. According to the thresholds specified in the text, most points fall into the mixed ER + CR state. The cluster at the excitonic resonance corresponds to the initial evolution (<200 fs).

Let us now examine the time evolution of the excitation as it progresses along the dynamics. At time zero, the Ω-matrix is dominated by large diagonal elements (Ω₁₁ and Ω₂₂), as shown in Fig. 11-top. Concurrently, the off-diagonal terms (Ω₁₂ and Ω₂₁) are nearly zero, consistent with a low-charge-transfer state. It implies that initially, [3.3]PCP is excited into an excitonic resonance state across both rings of the type .¹² To support this observation, the trajectories were each individually classified into the different types of states at each time step to avoid artifacts of the averaging of the Ω-matrix elements. The traces of fractions of trajectories evolving in each type of state are shown in Fig. 11-bottom. This corroborates the observation from the Ω-matrix elements, clearly showing that around 90% of the initial conditions correspond to an ER state. Around 10% are immediately excited to a mixed ER + CR state.

Fig. 11 (top) Ω elements of the fragment-based charge transfer matrix averaged over all trajectories as a function of time for the current (occupied) state. (bottom) Time evolution of populations of different excited state characters based on the classification of the individual trajectory at each time step.

The natural transition orbitals (NTO) of the third excited state (S₃) at the S₀ minimum geometry, where dynamics begin, are shown in Fig. 12a. The top three NTO pairs account for 34%, 31%, and 23%, respectively, totaling 88% of the total transition. The 1TDM descriptors for this structure closely match the mean values of the initial point in the dynamics: the CT value is ≈0.113, and the PR is 2.22.

Fig. 12 Main NTOs of [3.3]PCP for (a) S₃ state at the S₀ minimum (B_u-like in approximate C_2h symmetry) and (b) S₁ state (B_g-like) at S₁ minimum. Each case illustrates the hole, the electron, and their respective contributions (%) to the excited state.

Within the first 200 fs, a notable shift occurs: Ω₁₁ and Ω₂₂ decrease sharply to values around 0.22–0.25, while Ω₁₂ and Ω₂₁ increase significantly, reaching values around 0.15–0.18. This redistribution of Ω elements indicates the onset of inter-fragment electronic coupling, resulting in a transition from excitonic resonance to a state with increased charge-transfer character. The symmetry of the Ω-matrix elements, namely Ω₁₁ ≈ Ω₂₂ and Ω₁₂ ≈ Ω₂₁, means that the molecule is a state with mixed excitonic resonance and charge-resonance characters. This is supported through the individual classification of the trajectories. After 200 fs, almost 100% of the trajectories in an excited state occupy a mixed ER + CR state. Only tiny fluctuations can be observed at the end of the propagation time.

The character of the S₁ state at its minimum geometry was also analyzed using natural transition orbitals (NTOs), as shown in Fig. 12b. This analysis revealed a single dominant excitation, with a leading pair contributing 86% of the transition density. The values for CT (∼0.45) and PR (∼2.33) for this structure are nearly the same as those recorded during the dynamics from 200 fs to 2000 fs.

4. Conclusions

This study provides a detailed picture of the excited-state relaxation dynamics of [3.3]PCP, utilizing TDDFT and surface-hopping simulations in conjunction with electronic density analysis. Vertical excitations identified S₃ as the lowest-energy state with non-zero oscillator strength, which funnels within ∼0.1 ps via S₂ to S₁. Once in S₁, the molecule becomes kinetically trapped, with more than 85% of the trajectories remaining in this state throughout the 2-ps simulation window.

Taken together, the results outline the following photophysical mechanism. Excitation into an excitonic-resonant S₃ state is followed by ultrafast internal conversion to S₁, where the system stabilizes as a mixed excitonic-resonance/charge-resonance state delocalized over the two rings. In this regime, the dynamics are dominated by the inter-ring breathing motion—a low-frequency vibration also reported for [2.2]PCP—showing that both systems relax along a common excimer-forming coordinate.

Structurally, the dimer adopts a boat-type puckering at the para carbons, with these carbons bent toward the opposite ring. This geometry facilitates an approach to the S₁/S₀ crossing, but the S₁ → S₀ decay remains kinetically hindered, with an extrapolated lifetime of 12–23 ps. In condensed phases—solution or crystal—dissipation of excess kinetic energy would likely prolong ground-state recovery even further.

The combination of restricted structural distortion and sustained electronic delocalization explains the persistence of the excited state and underscores the ability of [3.3]PCP to sustain long-lived, delocalized excitations. At the same time, the complex S₁/S₂ coupling and slow S₁ → S₀ relaxation confirm that [3.3]PCP is a stringent test case for nonadiabatic dynamics schemes. While the surface hopping framework reproduces the key features of population transfer and trapping, a complete description of coherence and energy dissipation will require more advanced approaches. These characteristics make PCP scaffolds not only relevant for optoelectronic applications but also ideal for benchmarking theoretical methods for excited-state dynamics.

In this context, TDDFT-based surface hopping should be regarded as a approximate baseline: it captures the fast relaxation pathways and the emergence of long-lived S₁ trapping in [3.3]PCP, yet its quantitative predictions for the long S₁ lifetimes and S₁/S₀ decay remain method-dependent, especially in regions where multireference effects and true conical intersections become important, and are additionally affected by the use of the TD-BA approximation for nonadiabatic couplings in a dense low-lying manifold. The present work, therefore, establishes [3.3]PCP as a demanding test case that future nonadiabatic dynamics studies at higher levels of electronic-structure theory can build upon.

Author contributions

Conceptualization: MB; data curation: MTH; formal analysis: MTH, LMI; funding acquisition: MB; investigation: MTH; methodology: MTH, RSM, LMI; project administration: MB; software: RSM; supervision: MB, LMI; visualization: MTH, LMI; writing – original draft: MTH, LMI, MB; writing – review & editing: MTH, MB, LMI.

Conflicts of interest

There are no conflicts to declare.

Data availability

Simulated IR spectrum, additional surface-hopping results, geometry illustrations, the internal-conversion energy pathway, and Cartesian coordinates of the minima and the intersection point have been included in the supplementary information (SI). See DOI: https://doi.org/10.1039/d5cp04083c. Dynamics results are available at https://doi.org/10.5281/zenodo.15847715.

Acknowledgements

This work received support from the French government under the France 2030 investment plan as part of the Initiative d’Excellence d’Aix-Marseille Université (A*MIDEX AMX-22-IN1-48). The authors acknowledge the Centre de Calcul Intensif d’Aix-Marseille for granting access to its high-performance computing resources. Erasmus Mundus Joint Master Fellowships in Chemical Nanoengineering supported MTH.

References

1. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Chem. Rev., 2007, 107, 926–952,  DOI:10.1021/cr050140x.
2. J. Mei, N. L. Leung, R. T. Kwok, J. W. Lam and B. Z. Tang, Chem. Rev., 2015, 115, 11718–11940,  DOI:10.1021/acs.chemrev.5b00263.
3. H. Hopf, Angew. Chem., Int. Ed., 2008, 47, 9808–9812,  DOI:10.1002/anie.200800969.
4. M. Auffray, D. H. Kim, J. U. Kim, F. Bencheikh, D. Kreher, Q. Zhang, A. D'Aléo, J.-C. Ribierre, F. Mathevet and C. Adachi, Chem. – Asian J., 2019, 14, 1921–1925,  DOI:10.1002/asia.201900401.
5. H. Kaji, H. Suzuki, T. Fukushima, K. Shizu, K. Suzuki, S. Kubo, T. Komino, H. Oiwa, F. Suzuki, A. Wakamiya, Y. Murata and C. Adachi, Nat. Commun., 2015, 6, 8476,  DOI:10.1038/ncomms9476.
6. M. Y. Wong and E. Zysman-Colman, Adv. Mater., 2017, 29, 1605444,  DOI:10.1002/adma.201605444.
7. Q. Zhang, J. Li, K. Shizu, S. Huang, S. Hirata, H. Miyazaki and C. Adachi, J. Am. Chem. Soc., 2012, 134, 14706–14709,  DOI:10.1021/ja306538w.
8. F. Würthner, Acc. Chem. Res., 2016, 49, 868–876,  DOI:10.1021/acs.accounts.6b00042.
9. E. Fauvel, T. Moussounda Moussounda Koumba, F. El Kadiry, S. Maria, M. Rollet, M. Maresca, D. Siri, J.-L. Clément, D. Gigmes and M. Nechab, Angew. Chem., Int. Ed., 2025, 64, e202422253,  DOI:10.1002/anie.202422253.
10. H. Vazquez, R. Skouta, S. Schneebeli, M. Kamenetska, R. Breslow, L. Venkataraman and M. S. Hybertsen, Nat. Nanotechnol., 2012, 7, 663–667,  DOI:10.1038/nnano.2012.147.
11. O. Haggag, R. Baer, S. Ruhman and A. I. Krylov, J. Chem. Phys., 2024, 160, 124111,  DOI:10.1063/5.0196641.
12. F. Plasser and H. Lischka, J. Chem. Theory Comput., 2012, 8, 2777–2789,  DOI:10.1021/ct300307c.
13. G. Miao and J. Subotnik, J. Phys. Chem. A, 2019, 123, 5428–5435,  DOI:10.1021/acs.jpca.9b03188.
14. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57,  DOI:10.1016/j.cplett.2004.06.011.
15. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106,  DOI:10.1063/1.2834918.
16. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, J. Chem. Phys., 2020, 152, 184107,  DOI:10.1063/5.0004635.
17. F. Weigend, F. Furche and R. Ahlrichs, J. Chem. Phys., 2003, 119, 12753–12762,  DOI:10.1063/1.1627293.
18. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112,  DOI:10.1063/1.4993215.
19. I. Lyskov, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2016, 144, 034104,  DOI:10.1063/1.4940036.
20. S. Grimme and M. Waletzke, J. Chem. Phys., 1999, 111, 5645–5655,  DOI:10.1063/1.479866.
21. C. Hättig, Adv. Quantum Chem., 2005, 50, 37–60,  DOI:10.1016/S0065-3276(05)50003-0.
22. S. Battaglia and R. Lindh, J. Chem. Phys., 2021, 154, 034102,  DOI:10.1063/5.0030944.
23. A. Heil, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2018, 149, 164106,  DOI:10.1063/1.5050476.
24. G. Li Manni, I. Fdez Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, S. I. Bokarev, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. De Graaf, M. G. Delcey, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov, L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krośnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, P. López Ríos, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Mörchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. S. Oakley, M. Olivucci, M. Oppel, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, I. Rivalta, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, T. R. Scott, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, Y. Shu, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023, 19, 6933–6991,  DOI:10.1021/acs.jctc.3c00182.
25. B. G. Levine, J. D. Coe and T. J. Martínez, J. Phys. Chem. B, 2008, 112, 405–413,  DOI:10.1021/jp0761618.
26. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627,  DOI:10.1021/j100096a001.
27. R. Crespo-Otero and M. Barbatti, Theor. Chem. Acc., 2012, 131, 1237,  DOI:10.1007/s00214-012-1237-4.
28. M. Barbatti, M. Bondanza, R. Crespo-Otero, B. Demoulin, P. O. Dral, G. Granucci, F. Kossoski, H. Lischka, B. Mennucci, S. Mukherjee, M. Pederzoli, M. Persico, M. Pinheiro Jr, J. Pittner, F. Plasser, E. Sangiogo Gil and L. Stojanovic, J. Chem. Theory Comput., 2022, 18, 6851–6865,  DOI:10.1021/acs.jctc.2c00804.
29. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071,  DOI:10.1063/1.459170.
30. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114,  DOI:10.1063/1.2715585.
31. M. T. do Casal, J. M. Toldo, M. Pinheiro Jr and M. Barbatti, Open Res. Europe, 2021, 1, 49,  DOI:10.12688/openreseurope.13624.1.
32. K. K. Baeck and H. An, J. Chem. Phys., 2017, 146, 064107,  DOI:10.1063/1.4975323.
33. J. M. Toldo, R. S. Mattos, M. Pinheiro Jr, S. Mukherjee and M. Barbatti, J. Chem. Theory Comput., 2024, 20, 614–624,  DOI:10.1021/acs.jctc.3c01159.
34. D. Cremer and J. A. Pople, J. Am. Chem. Soc., 1975, 97, 1354–1358,  DOI:10.1021/ja00839a011.
35. J. C. A. Boeyens, J. Chem. Crystallogr., 1978, 8, 317–320,  DOI:10.1007/BF01200485.
36. C. A. G. Haasnoot, J. Am. Chem. Soc., 1992, 114, 882–887,  DOI:10.1021/ja00029a013.
37. F. Plasser, J. Chem. Phys., 2020, 152, 084108,  DOI:10.1063/1.5143076.
38. M. Pinheiro, M. de Oliveira Bispo, R. S. Mattos, M. Telles do Casal, B. Chandra Garain, J. M. Toldo, S. Mukherjee and M. Barbatti, Digital Discovery, 2025, 4, 666–682,  10.1039/D4DD00374H.
39. G. Braun, I. Borges Jr, A. Aquino, H. Lischka, F. Plasser, S. A. do Monte, E. Ventura, S. Mukherjee and M. Barbatti, J. Chem. Phys., 2022, 157, 154305,  DOI:10.1063/5.0113908.
40. S. Matsika and D. R. Yarkony, J. Am. Chem. Soc., 2003, 125, 10672–10676,  DOI:10.1021/ja036201v.
41. S. Mukherjee, R. S. Mattos, J. M. Toldo, H. Lischka and M. Barbatti, J. Chem. Phys., 2024, 160, 154306,  DOI:10.1063/5.0203636.
42. E. G. F. de Miranda, R. Souza Mattos, S. Mukherjee, J. M. Toldo, C. H. Choi, M. T. D. N. Varella and M. Barbatti, J. Chem. Theory Comput., 2025, 22, 1–19,  DOI:10.1021/acs.jctc.5c01529.
43. O. S. Haggag, N. Levinsky and S. Ruhman, ChemPhotoChem, 2022, 6, e202200181,  DOI:10.1002/cptc.202200181.
44. L. M. Ibele, P. A. Sánchez-Murcia, S. Mai, J. J. Nogueira and L. González, J. Phys. Chem. Lett., 2020, 11, 7483–7488,  DOI:10.1021/acs.jpclett.0c02193.

---