Elucidating the structure and binding nature of thianaphthene dimers using gas-phase infrared spectroscopy

Abstract

Crucial for many biological systems, in astrochemistry, and for fundamental chemistry in general, the conformations adopted by weakly bound complexes of polycyclic aromatic hydrocarbons (PAHs) have been the focus of debate for decades. Still, there are great challenges in forming such complexes in the laboratory, measuring accurate spectroscopic information sensitive to structure, and computing molecular properties for systems subtly bound by dispersive interactions. Here, we employ a combination of gas-phase infrared spectroscopy with extensive density functional theory (DFT) calculations to unambiguously determine the preferred conformation adopted by dimers of neutral thianaphthene, a PAH composed of one six- and one five-membered ring with an incorporated sulfur atom. A very wide spectral range from 350 to 3150 cm⁻¹ is covered, allowing a determination of the effect of complexation on fingerprint vibrations as well as C–H stretches. A comparison of the recorded infrared spectra of monomers and dimers in combination with detailed vibrational calculations assigns a π-stacked configuration for the complex. This agrees with energetic arguments, where DFT predicts the isomeric T-shaped configuration to be 0.13 eV higher in energy. The potential energy surface of the complex is explored using the nudged elastic band (NEB) method and the nature of the interaction between neutral monomers is investigated based on the local energy decomposition (LED) analysis. The π-stacked dimer is overwhelmingly stabilized by π⋯π dispersion, an interaction that is much weaker in the T-shaped configuration, despite the effect of C–H⋯π forces. The methodology applied here to thianaphthene is extendible to dimers without a permanent dipole moment, hence invisible to microwave spectroscopy, as well as larger clusters of PAHs.

---

Introduction

Polycyclic aromatic hydrocarbons (PAHs), a class of organic species composed of two or more aromatic rings fused together, have been the centre of research for decades. In an astrochemical context, PAHs are crucial, as this family of molecules is observed throughout a wide range of astronomical environments.¹ For example, PAHs are thought to be responsible for the so-called aromatic infrared bands (AIBs), a series of emission lines most prominently centred at 3.3, 6.2, 7.7, 8.6, 11.3, and 12.7 µm, common frequencies of the vibrational modes of these molecules.^2,3 Moreover, many substituted PAHs that possess a permanent dipole moment have been detected in space via their rotational spectra based on a combination of microwave spectroscopy and radioastronomy.^4–6 On Earth, PAHs are relevant for combustion chemistry, as they are crucial for soot formation and are considered environmental pollutants.⁷

From a chemical point of view, the non-covalent interaction in complexes of PAHs is of great significance, given their crucial role in many biological and chemical systems,^8,9 as well as providing unique data sets for benchmarking quantum chemical calculations.¹⁰ In this context, interest stems from the subtle nature of the interaction between closed-shell aromatic species, governed by an interplay of electrostatic forces, Pauli repulsion, and dispersion forces.¹¹ In PAHs clusters in general, but mostly investigated on dimers in particular, such interplay of forces leads primarily to two structural families, the so-called π-stacked and T-shaped configurations.¹² In the former, monomers bind primarily by a π⋯π interaction, leading to parallel conformations, whereas in the latter C–H⋯π forces are significant and dimers orient perpendicular to each other.¹³ In addition, substituted PAHs also support hydrogen bonding structures, with monomers oriented on the same plane.¹⁴ The benzene dimer, a prototypical system for studying the structure and non-covalent interactions on aromatic species, has been studied in detail, with conflicting results from quantum chemical calculations,^15–17 partly resolved by microwave spectroscopy measurements revealing a T-shaped configuration.¹⁸ However, the π-stacked structure of benzene does not possess a permanent dipole moment, making it undetectable with microwave spectroscopy. In addition, infrared spectroscopy, which does not suffer from that limitation, has been employed to investigate the conformation of dimers of pyrazine,¹⁹ pyridine,²⁰ azaphenanthrene,¹⁴ and naphthalene and acenaphthene,²¹ revealing system-specific conformations. Nevertheless, despite the work currently done, there is still a lack of widespread experimental information about dimers of aromatic species.

More recently, microwave spectroscopy was used to identify the conformation adopted by dimers of thiophene, an analogous species to benzene with a sulfur atom incorporated in its ring.²² In this case, thiophene does possess a permanent dipole moment, so the measurements are sensitive to both π-stacked and T-shaped structures. Nevertheless, similar to the early experiments on benzene, only a T-shaped configuration was identified. This raises the question about the effect of size on the preferred conformation of dimers of sulfur-containing PAHs (S-PAHs), since naphthalene, for example, adopts a π-stacked geometry,²³ thus different from the smaller benzene. This suggests a conformational change with PAH size. In S-PAHs, also the S⋯π interaction plays a role in determining the conformation of neutral dimers, an important non-covalent interaction in biological systems.²⁴ Moreover, S-PAHs are gaining attention for their possible presence in the interstellar medium, where they have been proposed as a possible sink of interstellar atomic sulfur,²⁵ thus contributing to the missing sulfur in the so-called sulfur depletion problem.²⁶ Notably, the sulfur-bearing PAH 2,5-cyclohexadien-1-thione has been recently discovered in space.²⁷

In this work, we employ a combination of gas-phase infrared spectroscopy and detailed density functional theory (DFT) calculations to identify the preferred configuration of thianaphthene dimers, an indene analogue with the CH₂ group replaced with the sulfur atom, thus the next size step from thiophene. The use of the FELIX free electron laser²⁸ allows the measurement of infrared spectra in a wide spectral range, covering the fingerprint region of PAHs, in addition to the C–H stretching range sampled with a pulsed infrared OPO/OPA laser. Our analysis conclusively shows that thianaphthene dimers adopt a π-stacked geometry, and provide a detailed understanding of the non-covalent interactions at play, leading to this preferred conformation.

Methods

Experimental

Experiments are conducted in a molecular beam setup coupled to a beamline of the FELIX free electron laser of the HFML-FELIX Institute (Nijmegen, The Netherlands). Briefly, a reservoir containing thianaphthene powder is resistively heated to 80 °C. The created vapor is then expanded supersonically into vacuum by a stream of Ar gas at a backing pressure of 2 bar via a pulsed valve, forming a molecular beam of neutrals. Following, the beam is collimated by a 2 mm skimmer before entering a perpendicularly extracted reflectron time-of-flight mass spectrometer, where the composition of the molecular beam is analysed after ionization. At the low vibrational temperatures of the molecular beam, estimated to be close to 40 K,²⁹ a wide distribution of thianaphthene clusters is formed, (Thia)_N with N spanning up to 9 on some cases. Here, however, experimental conditions are tuned to primarily produce the thianaphthene monomer and its dimer. To record mass spectra, the molecular beam is ionized via resonance-enhanced multiple photon ionization (REMPI) at the mass spectrometer entrance using the light of a tunable dye laser operating on Rhodamine B in ethanol, pumped by the second-harmonic of a Nd:YAG laser (532 nm; Spectra-Physics). After frequency doubling, the dye laser (Lioptec) is tuned in the 293–298 nm range, where the S₁ ← S₀ transition of thianaphthene has been observed.³⁰

To measure the infrared spectrum of (Thia)₁ and (Thia)₂, the counterpropagating light of FELIX is used to resonantly excite the gas-phase species vibrationally. This excitation depopulates their vibrational ground-state, making ionization less efficient and thus inducing a signal depletion in mass spectra (the so-called ion-dip spectroscopy).³¹ In order to account for signal fluctuations, FELIX runs at half the repetition rate of the molecular beam experiment, 10 and 20 Hz, respectively. This allows the recording of consecutive mass spectra with and without IR light. The IR yield is defined as I_IR = −ln(I/I₀)/P, with I and I₀ the intensity in mass spectra with and without FELIX excitation, respectively, and P the laser energy. FELIX is scanned in the so-called fingerprint region, ranging from 350 to 2200 cm⁻¹, in steps of 1 cm⁻¹. Importantly, FELIX is timed to interact with the molecular beam 300 µs prior ionization, thus ensuring that the recorded IR spectra correspond to species in the neutral state. The linewidth of FELIX is wavelength dependent, roughly corresponding to a FWHM of 0.5% of the central wavenumber in the current measurements. Further experimental details can be found elsewhere.^32,33

In addition to the fingerprint region, the range where the C–H stretching modes are located (3000–3150 cm⁻¹) is covered using a tabletop infrared OPO/OPA (LaserVision), with a higher spectral resolution of 0.1 cm⁻¹. Here, the IR radiation is crossed perpendicular to the molecular beam, opposite in direction to the ionization laser, and is timed 2 µs prior ionization. The OPO is scanned in steps of 0.16 cm⁻¹. Further details of this procedure can also be found elsewhere.³⁴

Computational

Density functional theory calculations (DFT) are conducted with the ORCA 6.01 software package.³⁵ For (Thia)₂, an initial structural search was performed using the GOAT algorithm,³⁶ preserving the topology of the thianaphthene monomers but varying their relative orientation in the dimer. This process yielded five distinct conformations, all corresponding to π-stacked geometries with different angles between the monomer units. Following, the set of structures are re-optimized at the DFT level, employing the meta-GGA TPSS functional³⁷ in conjunction with the Def2-TZVP basis set³⁸ and D3 dispersion corrections.³⁹ All π-stacked structures are found to be essentially isoenergetic, with relative energies below 0.001 eV. These structures are presented in the supporting information (SI). Furthermore, new configurations were manually constructed, in particular focusing on T-shaped geometries as well as side-by-side structures, where the two monomers bind horizontally forming C–H⋯C–H or C–H⋯S bonds. This process yielded one stable T-shaped geometry and no horizontal configuration. The DFT calculations are performed with the “tight” convergency criteria for the SCF and geometry optimization cycles, as implemented in ORCA, as well as with the finer “defgrid3” integration grid. For all stable structures, harmonic frequency calculations are carried out.

Anharmonic frequency calculations for the dimers are attempted but failed due to issues with the dispersion forces and large amplitude motions in computing the force constants on the potential energy and dipole moment surfaces. Nevertheless, for an accurate description of the vibrational spectrum of the thianaphthene monomer, anharmonic vibrational frequencies are computed using the B3LYP hybrid DFT functional⁴⁰ in conjunction with the aug-cc-pV(T+d)Z basis set.^41,42 First, the geometry of the monomer is optimized to the ground state minimum and the normal modes and harmonic frequencies are computed. Following this, a quartic force field (QFF) is computed at the same level of theory. A QFF is a fourth-order Taylor expansion of the potential energy surface surrounding the equilibrium geometry and consists of the quadratic, cubic, and semi-diagonal quartic force constants.⁴³ Anharmonic vibrational frequencies are then computed via second order vibrational perturbation theory (VPT2)^44,45 using Gaussian 16⁴⁶ including up to 3-quanta vibrational modes.

Moreover, of the different π-stacked configurations of the dimer, the lowest in energy is selected for further analysis, together with the T-shaped geometry. The lowest-energy path connecting these structures within the PES of the dimer is explored using the Nudged Elastic Band (NEB) method implemented in ORCA. For this, the same level of theory as for the DFT calculations is employed to compute 8 images along the energy path. The same methodology is used to compute the path along two π-stacked geometries. In addition, the interaction strength between monomers in both geometries is analysed using the local energy decomposition (LED) methodology implemented in ORCA, for which single-point DLPNO-CCSD(T) calculations are conducted using the structures optimized with DFT. This allows decomposing the interaction energy of the monomers into dispersive and non-dispersive terms, as discussed below and in the literature.^47,48

Results and discussion

Size distribution and REMPI spectra

The cold conditions of the molecular beam allow the formation of a distribution of thianaphthene clusters, (Thia)_N, as shown by the representative mass spectrum in Fig. 1(a). An identification of the true size-to-size abundances of the different (Thia)_N complexes in the molecular beam is hard to determine, given the strong sensitivity to ionization efficiency when recording mass spectra and the possibility of some fragmentation occurring upon ionization, but as shown in the figure, complexes with up to N = 6 are clearly present, while an extended mass spectrum shows even larger clusters (see the SI). Here, we focus on the monomer (Thia)₁ and dimer (Thia)₂ species, although it is noteworthy that such large complexes are stable in the molecular beam.

Fig. 1 (a) Mass spectrum of neutral thianaphthene clusters in the molecular beam, constructed as the average over all scanned REMPI resonances. (b) Typical REMPI spectra of the thianaphthene monomer, dimer and trimer, highlighting the wavelength used to measure the IR spectrum of the dimer.

Important for conducting spectroscopic measurements on such wide mass distributions is to have distinct REMPI resonances per complex, such that specific species can be isolated in the recorded mass spectra. Fig. 1(b) shows that this is indeed the case here, to a large extent. The S₁ ← S₀ transition of the monomer, seen at 293.52 nm, is sharp and at a slightly higher energy than the resonances of the dimer and trimer. For these species, the REMPI spectra are broad, as typically seen in weakly bound complexes,^14,29 and while they do overlap, selecting 296 nm as the excitation wavelength allows mostly the selective ionization of the dimer. Thus, the experiments discussed next were performed at 293.52 nm for the monomer and at 296 nm for the dimer. We note that in a previous study the S₁ ← S₀ transition of thianaphthene was reported as 293.58 nm, very close to our measurement.³⁰

Infrared spectroscopy of thianaphthene monomer

The measured IR spectrum for the yet-cold thianaphthene monomer is depicted in Fig. 2, showing in panel (a) the fingerprint region between 350 and 2300 cm⁻¹ (FELIX) and in panel (b) the C–H stretching range, between 3000 and 3150 cm⁻¹ (OPO). Panel (a) reveals a dense infrared spectrum, where the presence of sulfur breaks the molecular symmetry leading to all 38 vibrational modes being infrared active. The measured spectrum (in filled blue) agrees well with data from the NIST data base,⁴⁹ although here recorded with higher resolution due to the low temperature. Among the many features, the three most intense bands are found at 689, 732 and 761 cm⁻¹, with similar intensities, as well as other significant bands seen at 409, 467, 555, 799 and 867 cm⁻¹. In addition to the experimental spectrum, panel (a) presents in red the computed anharmonic vibrational transitions of thianaphthene at the B3LYP/aug-cc-pV(T+d)Z level of theory, revealing a very good match for most features, despite some small frequency differences and mismatches in relative intensities. We note that between 1600 and 2100 cm⁻¹ all observed bands are of anharmonic origin (combination bands and overtones, both 2- and 3-quanta transitions), a trend seen in many other PAH systems, both unsubstituted and substituted.^50,51 This explains why previous computational analysis of thianaphthene based on harmonic calculations did not reproduce that part of the IR spectrum.²⁵

Fig. 2 Infrared spectrum of thianaphthene, measured with FELIX in the 350–1700 cm⁻¹ range (a) and with the OPO/OPA laser in the 3000–3150 cm⁻¹ region (b). The red stick spectrum corresponds to the anharmonic vibrational transitions computed using the B3LYP/N07D method. The geometry of thianaphthene is shown as an inset in panel (a).

The three most intense features observed in Fig. 2a correspond to out-of-plane C–H bending modes. The IR spectrum of indene⁵² exhibits only two intense features in this frequency region at 717.8 and 765.7 cm⁻¹. The substitution of a sulfur atom for the sp³ hybridized CH₂ group leads to changes in the vibrational spectrum, most notably the presence of a third intense feature in thianaphthene. In indene, the strong transition at 718 cm⁻¹ corresponds to an asymmetric out-of-plane bending motion involving the four hydrogens on the 6-membered ring (ν₃₆) while the transition at 766 cm⁻¹ originates from a symmetric out-of-plane CH bending of the four CH groups on the 6-membered ring in addition to the CH group on the alpha carbon (ν₃₄). In thianaphthene, the asymmetric out-of-plane CH bending of the four 6-membered ring hydrogens (ν₂₈) occurs at 732 cm⁻¹, a blue shift of 13 cm⁻¹ from indene. Interestingly, the symmetric out-of-plane CH bending (761 cm⁻¹; ν₂₇) does not seem to be as affected by the sulfur substitution, with a redshift of only 5 cm⁻¹ in thianaphthene. The thianaphthene band at 689 cm⁻¹ is unique to this molecule and arises from anharmonic mixing. This feature consists of an out-of-plane CH bending of the two hydrogens on the 5-membered ring (ν₃₀) mixed with a combination band of a rocking motion and an asymmetric in-plane carbon skeleton deformation (ν₃₄ + ν₃₈). Additionally, several vibrational modes involving displacements of the C–S bond are identified in the IR spectrum, in particular those observed at 706, 796, 867, 1055 and 1093 cm⁻¹.

In addition to the mid-infrared range, Fig. 2b shows the C–H stretching region, covered with a tabletop infrared OPO/OPA. Given the low temperature of our gas-phase measurement, its resolution here is far greater than that of the NIST data,⁴⁹ allowing a clear distinction of individual IR features, in particular with intense bands at 3067, 3073 and 3086 cm⁻¹, together with many less intense modes. The stick spectrum in (b) shows the anharmonic vibrational frequencies, revealing a satisfactory agreement with the experiment. A well-known challenge for PAHs is computing C–H stretching modes, which often poorly reproduce experiments due to strong anharmonic effects.^53,54 Here, our calculations reproduce well the three main features of the experiment, as well as some of the weaker modes. Nevertheless, the match between experiment and theory is somewhat weaker than in the mid-infrared range, as expected.

Anharmonic vibrational frequency computations using various methodologies were performed on the thianaphthene monomer before arriving at the B3LYP/aug-cc-PV(T+d)Z level as the most accurate. The first attempt with the B3LYP/N07D level of theory, which has been shown to provide exceptionally accurate frequencies for standard and N-substituted PAHs,^{51,52,54–56} produced fairly good agreement in the mid-IR region, but performed quite poorly in the CH stretching region. The N07D basis set struggles for thianaphthene due to the underlying small 6-31G(d) double-ζ basis set. The inclusion of electron correlation using the rev-DSDPBEP86/jun-cc-pVDZ+B3LYP/N07D hybrid method⁵⁷ that accurately computes the vibrational frequencies of cyano-substituted PAHs^50,58 also struggled in the CH stretching region for thianaphthene. A switch to the Dunning family of basis sets and a size increase to the triple-ζ aug-cc-pV(T+d)Z basis set including the additional d functions is required to properly describe the presence of the 3rd period sulfur atom.

Dimer versus monomer infrared spectra

The presence of monomers and dimers within the same molecular beam allows the measurement of their infrared spectra under the same experimental conditions, meaning that IR spectra are sensitive to small variations upon complexation. This ultimately provides crucial information to elucidate the preferred conformation adopted by the dimer.^14,21 Fig. 3 compares the measured IR spectrum of the dimer to that of the monomer, with the former at the top and the latter at the bottom (multiplied by −1 for visualization purposes). Moreover, the figure is divided into panels (a) and (b), showing the mid-infrared and the C–H stretching regions at higher resolution, respectively.

Fig. 3 Comparison of the IR spectra of thianaphthene dimer and monomer, separated into the 350–1700 cm⁻¹ range in (a) and the 3000–3150 cm⁻¹ region in (b). For visualization purposes, the spectrum of the monomer is multiplied by −1. Vertical dashed lines are included to help visualize small differences in the spectra.

In the mid-infrared range, presented in panel (a), only small differences are visible between dimer and monomer. Close inspection, however, reveals important effects. For example, the relative intensity of the three most pronounced bands in the monomer change, with the feature at 761 cm⁻¹ (leftmost vertical dotted line) becoming roughly twice as strong in the dimer relative to the other bands in this region. Similarly, the relative intensity of the three peaks between 980 and 1100 cm⁻¹ is slightly affected by complexation. Moreover, some bands are redshifted in the dimer, as highlighted by the vertical dashed lines included in the figure (zoomed in versions are presented in the SI). For these two cases, for instance, redshifts of roughly 5 cm⁻¹ are observed in the dimer. Nevertheless, the IR spectrum of the dimer shares many similarities with that of the monomer, suggesting only a weak interaction in the complex.

Further differences are observed in the C–H stretching region, presented in Fig. 3(b). Here, redshifts of some bands are clear, in particular for the two strongest features. The redshifts are slightly smaller, close to 3 cm⁻¹, but significant if compared to the 0.1 cm⁻¹ resolution of the infrared OPO laser. Assuming a Gaussian shape, the experimental bands in this range have FWHM close to 1.5 cm⁻¹. Moreover, the spectrum of the dimer presents additional IR features in the spectrum, particularly around 3100 cm⁻¹, where the monomer only shows a sharp band at 3114 cm⁻¹, whereas for the dimer many peaks are seen between 3095 and 3115 cm⁻¹, partly overlapping. Similarly, between 3020 and 3050 cm⁻¹ three distinct peaks are observed for the monomer, in contrast to a much broader IR absorption for the dimer in that range. We stress that even though changes are subtle, measurements are conducted under the same excitation conditions, i.e., with the same FELIX settings, meaning that they cannot be attributed to experimental artifacts, such as different laser power during the measurements, or variations in laser wavelength calibration or linewidth.

Conformational analysis and structural assignment

In order to assign a conformation to (Thia)₂ based on the measured IR spectra, calculations of vibrational frequencies are required. As detailed in the introduction, two families of conformations are expected for the complex, namely the so-called π-stacked and T-shaped geometries. As discussed, our conformational search confirms this expectation, yielding 5 distinct π-stacked structures and a single T-shaped geometry. Within the π-stacked family, the relative orientation between monomers varies, although TPSS DFT calculations reveal quasi isoenergetic structures. Nevertheless, we select the lowest energy one for the analysis further on, while all structures are presented in the SI.

Fig. 4(a) depicts the most stable π-stacked geometry together with the single stable T-shaped conformation. At the TPSS DFT level, the T-shaped conformation lies 0.13 eV above the π-stacking geometry. Notwithstanding, TPSS DFT calculations on the smaller thiophene dimer give almost isoenergetic T-shaped and π-stacked conformations, with the T-shaped just lower by 0.001 eV, in agreement with previous studies.²² Thus, purely based on TPSS DFT energies, increasing the size of the S-PAH favours a π-stacking conformation. In Fig. 4(b), the minimum-energy path along the PES connecting both conformers is presented (blue circles), based on NEB calculations with ORCA. As seen, an energy barrier of 0.28 eV with respect to the lowest-in-energy π-stacked conformer (right side in the plot), is computed, a value too high to be accessible at 40 K, the typical vibrational temperature of the molecular beam. But even when assuming a Boltzmann distribution, an energy difference of 0.13 eV yields a population solely composed of the π-stacked conformer. In contrast, exploring the path between two π-stacked structures, as shown in Fig. 4(b) by red squares, gives the much lower energy barrier of 0.07 eV. This value is still high in comparison with 40 K, but within reach depending on the timescale. Therefore, the possibility of an ensemble of π-stacked conformers cannot be excluded. Moreover, the vibrational spectra of the π-stacked conformers are too similar as to distinguish them in the experiment, as shown in the SI.

Fig. 4 (a) Geometries of the thianaphthene dimer on the π-stacking and T-shaped conformations. Two representative distances between atoms in both structures are depicted. (b) PES along the minimum-energy path connecting the conformations shown in panel (a), presented by blue circles. The initial and final energies are depicted, as well as the located transition state. In addition, a PES scan along the path connecting two π-stacked conformations is presented by red squares.

Despite the strong computational conclusions based on energetics, an elucidation of the preferred conformation in the molecular beam requires experimental evidence. For this purpose, the harmonic vibrational frequencies of the geometries shown in Fig. 4(a) were computed, allowing a direct comparison with the measured IR spectrum of the thianaphthene dimer. The computed vibrational spectrum of the dimer on both conformations is presented in Fig. 5(a), covering the entire measured range from 350 to 3150 cm⁻¹, with the experimental data shown at the top of the panel. A careful inspection of the spectrum reveals key differences between conformers, as highlighted in panels (b) and (c) where the spectra of the two unique dimer configurations are compared with the spectrum of the monomer (also harmonic).

Fig. 5 (a) Computed vibrational spectra of the thianaphthene dimer in the π-stacked (green) and T-shaped (purple) conformations, with the experimental data presented at the top of the figure. Expanded views of the regions between 600–800 and 3050–3110 cm⁻¹ are presented in panels (b) and (c), respectively, which also present the spectrum of the monomer for comparison.

Fig. 5(b) focuses on the region between 650 to 800 cm⁻¹, corresponding to the C–H out-of-plane bending modes. As seen in comparison with the monomer, small redshifts are predicted for the π-stacked conformation, whereas a clear splitting and an overall blueshift of the three main modes are computed for the T-shaped geometry. The origin of the splitting is caused by the fact that for the T-shaped structure both monomers become distinguishable, with one pointing its C–H bonds towards the second, making the C–H out-of-plane bending modes of both units different. Clearly, this effect is not present for the π-stacked conformation, where the C–H out-of-plane bending modes of both monomers are only slightly destabilized. Thus, calculations are consistent with a π-stacked conformation for the dimer in the molecular beam, where the measured infrared spectrum presented in Fig. 3(a) revealed very similar spectra for the monomer and dimer, with only small redshifts of the main features upon complexation.

Similarly, the expanded view of the region of the C–H stretches (3050–3110 cm⁻¹), presented in Fig. 5(c), shows a clear splitting of the main modes in the T-shaped conformation, in addition to an overall blueshift of the features with respect to the monomer. In contrast, the vibrational spectrum of the π-stacked geometry is very similar to the monomer, only with small redshifts of the dominant peaks. Again, this aligns well with the experimental observations in this wavelength range (Fig. 3(b)), where the infrared spectrum of the dimer only showed small redshifts of some features, with an overall similar structure as the monomer spectrum. Therefore, our strategy of combining infrared spectroscopy measurements in molecular beams with vibrational frequency calculations based on DFT assigns with confidence a π-stacked conformation as the spectral carrier for the neutral thianaphthene dimer, stressing that the conformation can be a family of π-stacked geometries.

Intermolecular interactions and chemical bonding

Our combined experimental and computational analysis identifies a π-stacked configuration as the preferred orientation for the neutral thianaphthene dimer, in contrast to the smaller thiophene, shown by microwave spectroscopy to prefer a T-shaped geometry. This conclusion agrees with calculations of relative energies, which at the DFT level predict a more stable π-stacked conformation for thianaphthene (by 0.13 eV), whereas for the thiophene dimer both configurations have similar energies, with a slight preference for the T-shaped. In order to better understand the binding nature in the dimers, as well as to explore the change from preferred T-shaped to π-stacking conformation with S-PAH size, a local energy decomposition (LED) analysis is performed for the interaction energy in thianaphthene and thiophene dimers in the π-stacking and T-shaped conformations.

In the LED analysis, the total interaction energy is computed at the accurate DLPNO-CCSD(T) level of theory and is decomposed into different physically relevant components, as is explained in detail elsewhere.^47,48 In short, non-dispersive terms include: (i) the electronic preparation energy, corresponding to the energy required to perturb the wavefunction of the monomers from their ground-state configuration to the one they adopt in the dimer. (ii) The electrostatic interaction between monomers. (iii) The exchange interaction. (iv) The charge transfer energy. Of these four terms, the electronic preparation energy is always repulsive, while all others are attractive. In addition, the total energy includes a repulsive term denoted as strain, resulting from the distortion of the geometry of the monomers into their configuration adopted in the dimers. Finally, the total energy includes a stabilizing component accounting for the London dispersion attraction between monomers.^59,60 Fig. 6(a) summarizes the results from the LED analysis of (Thia)₂, where single-point DLPNO-CCSD(T) calculations are performed employing the geometry optimized with DFT.

Fig. 6 Local energy decomposition analysis of the π-stacked and T-shaped conformations of thianaphthene (a) and thiophene (b) dimers. For thianaphthene the geometries depicted in Fig. 4(a) are employed, whereas for thiophene the geometries presented in ref. 22 are used. (c) Dispersion interaction density plot of the thianaphthene dimer (isovalue = 0.12 kJ mol⁻¹ Bohr⁻³).

First, we note the negligible role played by strain, only accounting for 0.01 eV in the π-stacked dimer, whereas no strain is calculated for the T-shaped geometry. Therefore, the geometry of the monomers is slightly distorted, but only when the dimers adopt a π-stacked configuration; however, the effect is insignificant. Moreover, we see the dominance of dispersion, especially for the π-stacked configuration. Here, the dispersion term accounts for an attractive −0.41 eV, whereas the non-dispersion energy becomes repulsive, because the negative electrostatic, exchange and charge transfer terms are overcompensated by a large electronic preparation energy. Thus, the dimer is overwhelmingly stabilized by π⋯π dispersion. In contrast, although the dispersion term in the T-shaped geometry is still the largest, it is much lower than in the π-stacked configuration, only accounting for −0.22 eV, almost half the value in the π-stacked dimer. Moreover, because the repulsive electronic preparation is lower and the C–H⋯π electrostatic energy is larger, the non-dispersive term is also attractive. Hence, both dispersive and non-dispersive terms account for the interaction energy in the T-shaped geometry.

A visualization of the extent of the dispersion interaction in both configurations is found by plots of the dispersion interaction density (DID),⁶¹ as shown in Fig. 6(c). In the π-stacked conformation, dispersion is distributed along the monomers, with a strong involvement of the sulfur atom. In contrast, the T-shape geometry displays a weaker dispersion interaction, primarily due to the two C–H groups of the thianaphthene molecule pointing towards the second monomer. Thus, the DID plots agree with the LED analysis, showing a much stronger dispersion interaction in the π-stacked structure.

Ultimately, the total interaction energy is an interplay between the dispersive and non-dispersive terms. In the π-stacked dimer, dispersion is strong, but there is also an energy penalty for distorting the stable π clouds of thianaphthene. This effect is much weaker in the T-shaped configuration, given that only one of the monomers is affected. Thus, there is a subtle balance between the stabilizing π⋯π dispersion and the cost of distorting the stable electron density of the S-PAH. For thianaphthene, dispersion wins, and the dimer adopts a π-stacked conformation, as observed experimentally. Nevertheless, it is noteworthy that at the DLPNO-CCSD(T) level the preference for π-stacking over T-shaped is smaller than for DFT, only by 0.01 eV, showing the difficulties of this type of calculations.

Similarly, dispersion is the most significant component in the interaction energy of the thiophene dimer, as presented in Fig. 6(b), but its strength is weaker than for thianaphthene. Moreover, in stark contrast to thianaphthene, the difference in dispersion between π-stacked and T-shaped geometries is small, being only slightly larger for the former configuration. Still, the T-shaped dimer is stabilized by non-dispersive terms, whereas the energy penalty for electronic preparation in the π-stacked dimer compensates for the somewhat larger dispersion energy, thus favouring the T-shaped geometry. At the DFT level, the π-stacked and T-shaped geometries are almost isoenergetic, with just a small preference for the latter. Instead, at the DLPNO-CCSD(T) level the T-shaped is clearly preferred, being 0.06 eV lower in energy. Therefore, increasing the S-PAH size enhances the dispersion interaction between dimers, which overcomes the energy penalty for electronic preparation and triggers a transition for a T-shaped dimer in thiophene to a π-stacked one in thianaphthene.

Conclusions

The preferred conformation adopted by dimers of thianaphthene in the cold and isolated conditions of a molecular beam was unambiguously identified by combining infrared spectroscopy and density functional theory calculations. The neutral dimers adopt a π-stacked conformation, in contrast to the smaller thiophene dimer, which previous microwave spectroscopy experiments showed prefers a T-shaped configuration. DFT calculations support the assigned geometry of the thianaphthene dimer, predicting the π-stacked conformation 0.13 eV lower in energy than the T-shaped. A local energy decomposition analysis of the interaction energy of the dimer was performed at the accurate DLPNO-CCSD(T) level of theory, revealing a major influence of π⋯π dispersion in the π-stacked dimer, which more than compensates for a significant electronic preparation energy penalty. In the T-shaped geometry, in contrast, non-dispersive terms are attractive, but nonetheless, because of a much weaker dispersion interaction the conformer is less stable. Overall, dispersion is weaker in the smaller thiophene dimer, thus favouring a T-shaped conformation. Our work shows the accuracy of gas-phase infrared spectroscopy to infer the geometry of weakly bound complexes of PAHs, here applied to S-PAH, but nevertheless extendible to dimers without a permanent dipole moment, which therefore, cannot be investigated by microwave spectroscopy. Moreover, the same strategy can be employed to investigate larger clusters of PAHs.

Author contributions

G. Z. and P. F. performed the experiments. G. Z., P. F. and V. J. E. run the calculations. G. Z., P. F. V. J. E and S. B. wrote the manuscript. P. F. and S. B. envisioned the research.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are deposited in the Radboud Data Repository.

Supplementary information (SI): full mass spectrum of thianaphthene clusters, π-stacked conformations and relative energies, computed vibrational spectra of the π-stacked conformations, comparison of IR spectra of monomer and dimer. See DOI: https://doi.org/10.1039/d5cp01442e.

Acknowledgements

The authors gratefully acknowledge the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for the support of the HFML-FELIX Institute and for CPU time on the Dutch National Supercomputer Snellius (Project 2024.009). This work was possible thanks to the constant support of the FELIX staff. V. J. E. acknowledges Chapman University for financial support. This work used Bridges-2 at the Pittsburgh Supercomputing Center through allocation CHE250086 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

References

1. A. G. G. M. Tielens, Annu. Rev. Astron. Astrophys., 2008, 46, 289–337.
2. R. Chown, et al., Astron. Astrophys., 2024, 686, A75.
3. B. Khan, B. Abbott, E. Peeters, A. G. G. M. Tielens, T. Onaka, J. Cami, B. Schefter, C. Boersma, E. Dartois, J. R. Goicoechea and A. Maragkoudakis, Astron. Astrophys., 2025, 699, A133.
4. G. Wenzel, I. R. Cooke, P. Bryan Changala, E. A. Bergin, S. Zhang, A. M. Burkhardt, A. N. Byrne, S. B. Charnley, M. A. Cordiner, M. Duffy, Z. T. P. Fried, H. Gupta, M. S. Holdren, A. Lipnicky, R. A. Loomis, H. T. Shay, C. N. Shingledecker, M. A. Siebert, D. A. Stewart, R. H. J. Willis, C. Xue, A. J. Remijan, A. E. Wendlandt, M. C. McCarthy and B. A. McGuire, Science, 2024, 386, 810–813.
5. G. Wenzel, T. H. Speak, P. H. Changala, R. J. H. Willis, A. M. Burkhardt, S. Zhang, E. A. Bergin, A. N. Byrne, S. B. Charnley, Z. T. P. Fried, H. Gupta, E. Herbst, M. S. Holdren, A. Lipnicky, R. A. Loomis, C. N. Shingledecker, C. Xue, A. J. Remijan, A. E. Wendlandt, M. C. McCarthy, I. R. Cooke and B. A. McGuire, Nat. Astron., 2025, 9, 262–270.
6. G. Wenzel, S. Gong, C. Xue, P. Bryan Changala, M. S. Holdren, T. H. Speak, D. Archie Stewart, Z. T. P. Fried, R. H. J. Willis, E. A. Bergin, A. M. Burkhardt, A. N. Byrne, S. B. Charnley, A. Lipnicky, R. A. Loomis, C. N. Shingledecker, I. R. Cooke, M. C. McCarthy, A. J. Remijan, A. E. Wendlandt and B. A. McGuire, ApJL, 2025, 984, L36.
7. M. Thomson and T. A. Mitra, Science, 2018, 361, 978–979.
8. C. F. Matta, N. Castillo and R. J. Boyd, J. Phys. Chem. B, 2006, 110, 563–578.
9. R. Chelli, F. L. Gervasio, P. Procacci and V. Schettino, J. Am. Chem. Soc., 2002, 124, 6133–6143.
10. M. Rapacioli, F. Calvo, F. Spiegelman, C. Joblin and D. J. Wales, J. Phys. Chem. A, 2005, 109, 2487–2497.
11. S. Grimme, Angew. Chem. Int. Ed., 2008, 47, 3430–3434.
12. R. Podeszwa and K. Szalewicz, Phys. Chem. Chem. Phys., 2008, 10, 2735–2746.
13. R. Calinsky and Y. Levy, J. Phys. Chem. B, 2024, 128, 8687–8700.
14. X. Miao, T. Preitschopf, F. Sturm, T. Fischer, A. K. Lemmens, M. Limbacher and R. Mitric, J. Phys. Chem. Lett., 2022, 13, 8939–8944.
15. T. Sato, T. Tsuneda and K. Hirao, J. Chem. Phys., 2005, 123, 104307.
16. K. Carter-Fenk and J. M. Herbert, Chem. Sci., 2020, 11, 6758–6765.
17. C. R. Martinez and B. L. Iverson, Chem. Sci., 2012, 3, 2191–2201.
18. M. Schnell, U. Erlekam, P. R. Bunker, G. von Helden, J.-U. Grabow, G. Meijer and A. van der Avoird, Angew. Chem., Int. Ed., 2013, 52, 5180–5183.
19. J.-Y. Feng, Y.-P. Lee, P.-J. Hsu, J.-L. Kuo and T. Ebata, J. Phys. Chem. A, 2023, 127, 4291–4301.
20. O. Hübner, T. Thusek and H.-J. Himmel, Angew. Chem., Int. Ed., 2023, 62, e202218042.
21. A. K. Lemmens, P. Chopra, D. Garg, A. L. Steber, M. Schnell, W. J. Buma and A. M. Rijs, Mol. Phys., 2021, 119, e1811908.
22. X. Chen, G. Wang, X. Zeng, W. Li and M. Zhou, J. Am. Chem. Soc., 2025, 147, 28684–28693.
23. N. O. Dubinets, A. A. Safonov and A. A. Bagaturyants, J. Phys. Chem. A, 2016, 120, 2779–2782.
24. R. J. Zauhar, C. L. Colbert, R. S. Morgan and W. J. Welsh, Biopolymers, 2000, 53, 233–248.
25. X. J. Yang, L. Hua and A. Li, ApJ, 2024, 974, 30.
26. P. Ferrari, G. Berden, B. Redlich, L. B. F. M. Waters and J. M. Bakker, Nat. Commun., 2024, 15, 5928.
27. M. Araki, M. Sanz-Novo, C. P. Endres, P. Caselli, V. Rivilla, I. Jiménez-Serra, L. Colzi, S. Zeng, A. Megías, A. López-Gallifa, A. Martínez-Henares, D. San Andrés, S. Martín, M. A. Requena-Torres, J. García de la Concepción and V. A. Lattanzi, Nat. Astron., 2026, 10, 401–409.
28. D. Oepts, A. van der Meer and P. van Amersfoort, Infrared Phys. Technol., 1995, 36, 297–308.
29. A. K. Lemmens, P. Ferrari, D. Loru, G. Batra, A. L. Steber, B. Redlich, M. Schnell and B. Martinez-Haya, J. Phys. Chem. Lett., 2023, 14, 10794–10802.
30. M. R. Padhye and J. C. Patel, Trans. Faraday Soc., 1953, 49, 1119–1121.
31. S. Bakels, M.-P. Gaigeot and A. M. Rijs, Chem. Rev., 2020, 120, 3233–3260.
32. P. Ferrari, A. K. Lemmens and B. Redlich, Phys. Chem. Chem. Phys., 2024, 26, 12324–12330.
33. V. J. Esposito, P. Ferrari, W. J. Buma, R. C. Fortenberry, C. Boersma, A. Candian and A. G. G. M. Tielens, J. Chem. Phys., 2024, 160, 114312.
34. V. J. Esposito, P. Ferrari, W. J. Buma, C. Boersma, C. J. Mackie, A. Candian, R. C. Fortenberry and A. G. G. M. Tielens, Mol. Phys., 2024, 122, e2261570.
35. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
36. B. de Souza, Angew. Chem., Int. Ed., 2025, 64, e202500393.
37. V. K. Staroverov, J. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
38. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
39. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
40. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
41. J. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
42. J. T. H. Dunning, K. A. Peterson and A. K. Wilson, J. Chem. Phys., 2001, 114, 9244–9253.
43. R. C. Fortenberry and T. J. Lee, Computational vibrational spectroscopy for the detection of molecules in space, Annual Reports in Computational Chemistry, 2019, vol. 15, ch. 6, pp. 173–202.
44. J. K. G. Watson, Vibrational Spectra and Structure, Elsevier, Amsterdam, 1997.
45. I. Mills, Molecular Spectroscopy, Modern Research, Academic Press, New York, 1972.
46. M. J. Frisch, et al., Gaussian 16, Gaussian Inc., Wallingford CT, 2016.
47. G. Bistoni, WIREs Comput. Mol. Sci., 2020, 10, e1442.
48. A. Altun, R. Izsák and G. Bistoni, Int. J. Quantum Chem., 2021, 121, e26339.
49. https://webbook.nist.gov/cgi/cbook.cgi?ID=C95158&Mask=80#IR-Spec.
50. V. J. Esposito, P. Ferrari, C. Z. Palmer, C. Boersma, A. Candian, R. C. Fortenberry, W. J. Buma and A. G. G. M. Tielens, J. Phys. Chem. Lett., 2025, 16, 1296–1304.
51. C. J. Mackie, A. Candian, X. Huang, E. Maltseva, A. Petrignani, J. Oomens, A. L. Mattioda, W. J. Buma, T. J. Lee and A. G. G. M. Tielens, J. Chem. Phys., 2016, 145, 084313.
52. V. J. Esposito, S. Bejaoui, B. E. Billinghurst, C. Boersma, R. C. Fortenberry and F. Salama, MNRAS, 2024, 535, 3239–3251.
53. A. K. Lemmens, A. M. Rijs and W. J. Buma, ApJ, 2021, 923, 238.
54. E. Maltseva, A. Petrignani, A. Candian, C. J. Mackie, X. Huang, T. J. Lee, A. G. G. M. Tielens, J. Oomens and W. J. Buma, ApJ, 2015, 814, 23.
55. V. Barone, P. Cimino and E. Stendardo, J. Chem. Theory Comput., 2008, 4, 751–764.
56. V. Barone, M. Biczysko and J. Bloino, Phys. Chem. Chem. Phys., 2014, 16, 1759–1787.
57. R. Xu, Z. Jiang, Q. Yang, J. Bloino and M. Biczysko, J. Comput. Chem., 2024, 45, 1846–1869.
58. Y. Xu and M. Biczysko, Front. Chem., 2024, 12, 1439194.
59. I. Harden, F. Neese and G. Bistoni, Chem. Sci., 2023, 14, 10580.
60. G. Bistoni, A. Altun, Z. Wang and F. Neese, Acc. Chem. Res., 2024, 57, 1411–1420.
61. A. Wuttke and R. A. Mata, J. Comput. Chem., 2017, 38, 15–23.

---