Stacking structure in liquid polyaromatic hydrocarbons

Abstract

We present a neutron total scattering study of the polyaromatic hydrocarbons phenanthrene and pyrene, analysed by simulation based data refinement. We compare the liquid structure in terms of the degree and relative position of π-stacked pairs, to measurements of benzene and naphthalene. We find that pyrene has the highest level of parallel stacking, with 46% of molecules in stacked pairs. We observe a trend of decreasing displacement in the stacking with increasing aromatic core size, and find that for pyrene no offset is observed. The results have important implications for the fundamental understanding and modelling of polyaromatic hydrocarbon structure, and for applied fields including optoelectronics, astrochemistry and crude oil phase behaviour.

---

The non-covalent interactions between aromatic molecules and functional groups are fundamental to many biological, chemical and physical structures and processes.¹ Aromatic molecules with one or more fused rings are known as polyaromatic hydrocarbons (PAHs), for example naphthalene, phenanthrene and pyrene (Fig. 1). Their larger size tends to result in stronger intermolecular interactions, with a large flat area aligned for maximum dispersion attraction. The resulting structures are important for a wide range of science and technology, examples include understanding soot formation from combustion,² the aggregation of PAHs in the interstellar medium³ where the presence of clustering is important in interpreting the IR emission spectra,⁴ asphaltene aggregation in petroleum extraction and refining,⁵ and for organic electronics and optical devices where pyrene is a benchmark molecule, said to be the “fruit fly” for photochemists, with higher stacking resulting in higher charge-carrier mobilities.⁶ Results of ab initio calculations of pyrene dimers at SCC-DTFB,⁷ TPSS-D3BJ/def2-TZVPP⁸ levels of theory in vacuum and at ωB97X-D4/def2-TZVP level in toluene solution,⁹ show that offset or rotated parallel stacked dimers have the lowest energy, with only small differences in energy observed between different levels of offset, “graphite like” and rotated structures. Cluster calculations have been performed for pyrene using SCC-DTFB,⁷ showing that, for up to 3 or 4 molecules, a single stack is formed, with larger clusters being formed of multiple stacks. Phenanthrene and pyrene clusters have also been studied by MD simulation with the OPLS-AA force-field,¹⁰ with phenanthrene clusters taking herringbone structures, whereas pyrene clusters only showed amorphous structures.

Fig. 1 Chemical structures of benzene and PAH's with increasing numbers of aromatic rings and levels of condensation: naphthalene, phenanthrene and pyrene.

In this communication we report an experimental measurement of the liquid structures of phenanthrene and pyrene, which we combine with our previous analyses of liquid naphthalene¹¹ and benzene,¹² allowing us to study in detail changes in local structure with increasing number of aromatic rings and condensation. Measurement in the liquid state reduces the symmetry and packing constraints present in the crystal state, but still includes the complexities of multi-body interactions that are absent in dimer studies. It is therefore an excellent window into likely structural motifs in the liquid, solution and amorphous states where PAHs often have most applied relevance. In addition, these results provide benchmark experimental data for the basic understanding of the complex intermolecular interactions between aromatic molecules. In particular the nature of π stacked interactions is hotly debated^13–15 with subtle charge penetration effects playing a key role and the ultimate cause of horizontal displacements in stacked structures being ascribed to either dispersion attraction and exchange repulsion,¹⁴ or electrostatic effects.¹⁵

A detailed description of the experimental methods is presented at the SI. In summary, we measure the total neutron structure factor, F(Q) (equation S1), for 3 isotopologues of liquid phenanthrene (at 118 °C) and pyrene (at 160 °C): hydrogenated (H), deuterated (D) and a 1 : 1 molar mix of the two (HD). We obtain detailed information on the liquid structure using the simulation based refinement method, empirical structure refinement (EPSR),¹⁶ described in detail in SI Section S1. The refined simulated ensemble of structures is collected and analysed to give average structural information for the liquid.

Fits to the experimental F(Q) for the refined EPSR simulations are shown in Fig. 2 for pyrene and phenanthrene. The fit to the data is excellent, with some slight deviation for H containing samples, due to the well know difficulties in unambiguously removing inelastic scattering artifacts from samples containing light hydrogen. The clear improvement in fit from a simulation without refinement is shown in Fig. S3 and S4 in the SI. Fits to the real space total pair distribution function, f(r) (eqn (S3) and (S4)), are also excellent (SI Fig. S1 and S2). This function more clearly demonstrates the level of fit for short range correlations, therefore verifying the choice of bond lengths and energy parameters for the bonded potentials as outlined in the Methods section in the SI.

Fig. 2 Measured (dots) and EPSR simulated (red line) F(Q) for three isotopologues of liquid phenanthrene (a) and liquid pyrene (b).

We can now move forward with confidence that the refined simulations are consistent with the experimental data. Furthermore we note that EPSR is maximum entropy method, producing the most disordered, and therefore most likely structural ensemble that is not only constrained by the fit to the data, but also known quantities such as the molecular structure, the liquid density, and the calculated charge distribution across the molecule. Nevertheless, it is important to understand the possibility of alternative solutions to the fit. We can explore this by running data refinement with an alternative set of reference potentials using a “standard” set of fixed OPLS charges and Lennard-Jones parameters (as used for liquid naphthalene),¹¹ as described in Section S3 of the SI. Fig. S5 and S6 in the SI show how the quality of the fit is almost identical for the two different reference potential partial charge schemes, however some small structural differences are observed, as noted below.

The simplest structural parameter we obtain from the refined simulations is the one-dimensional radial distribution function, g(r), between molecular centres. We plot these for phenanthrene and pyrene, alongside those for liquid benzene (at 10 °C)¹² and naphthalene (at 85 °C),¹¹ in Fig. 3a). It is clear from the centre-centre g(r)'s that the local environment of phenanthrene and pyrene molecules are relatively unstructured (lower peak heights), but more complex than benzene, showing bumps and shoulders either side of the main peak in the g(r). Pyrene in particular shows a clear shoulder at r ≈ 4 Å. At this short distance, molecular geometry dictates that this must be parallel stacking of molecules. Phenanthrene is even less structured when looking at the g(r), likely due to the reduced symmetry of the molecule allowing a richer variety of nearest neighbour structures.

Fig. 3 (a) Radial distribution function g(r) between molecule centres for liquid aromatics phenanthrene and pyrene compared to benzene¹² and naphthalene¹¹ (b) Part of angular radial distribution function g(r,θ) for parallel stacked molecules only (θ = 0–10°) for the same aromatic liquids.

A key structural metric for aromatics is the angular radial distribution function, g(r,θ) (eqn (S5)), i.e. the g(r) as a function of the angle between the normal of the aromatic planes, θ. Our previous studies show a large increase in the probability of stacked interactions (low θ) when going from benzene¹² to naphthalene.¹¹ Full angular radial distribution functions are plotted in Fig. S7 and S8 in the SI, showing, in similarity to naphthalene, a clear preference for parallel stacked interactions for both phenanthrene and pyrene. In Fig. 3b we compare the stacked interactions between the different liquids, plotting g(r,θ = 0–10°) for each liquid on the same plot. Pyrene clearly shows the strongest preference for parallel stacked nearest neighbour interactions, with a strong peak at ∼3.7 Å (cf. graphite interlayer spacing of 3.35 Å). Pyrene also shows a second higher r peak at r = 7.8 Å, at roughly double the stacking peak distance, suggesting the presence of stacked trimers. This higher level of stacking in pyrene than smaller aromatics seems reasonable in the context of larger available overlap affording greater dispersion attraction. We note that this level of stacking for pyrene is somewhat affected by the choice of reference Coulomb potential, as shown in SI Fig. S10, the OPLS fixed charge potentials do show higher parallel stacking peak compared to the LigParGen calculated charges. The overall shape is still the same, but the peak at g(r,θ = 10°) is ∼40% lower. Therefore we caution that this region of g(r,θ) is most likely to be influenced by choices in the initial seed potential. For phenanthrene the parallel stacking peak at θ = 0–10° is considerably broader than for pyrene or even naphthalene, indicating much more disordered, and generally longer distance parallel interactions than seen in naphthalene and pyrene, likely due to the decreased symmetry of the molecule. For phenanthrene, no significant difference is observed in the g(r,θ) for the two partial charge schemes used (Fig. S9).

An interesting question now arises as to how many molecules in the liquid can be considered to be in stacked dimers. We define limits for r and θ from which to integrate g(r,θ), setting, θ = 0–20° and for benzene set an r_max = 5 Å due to the shoulder in g(r,θ = 0–10°) centred at ∼4 Å. From this we calculated an r_max for each liquid by scaling the benzene cut-off by the ratio of the cube-root of the molecular volume to account for differences in molecular density. Table 1 gives these cut-off's and the calculated stacking coordination numbers. The calculated cut-off's seem to be broadly justified, coming in close to, if slightly below the minima observed in the g(r,θ = 0–10°) for naphthalene and pyrene. Pyrene clearly has the largest number of stacked interactions, with close to half of all molecules being in a stacked pair. This falls to approximately a third for naphthalene and phenanthrene and and eighth of benzene molecules.

Table 1 Comparison of a stacking coordination number in each of the liquid aromatics as calculated by integrating g(r,θ) between θ = 0–20° and r = 0–r_max

	r max cut-off [Å]	Stacking coordination number
Benzene	5.0	0.126
Naphthalene	5.7	0.351
Phenanthrene	6.2	0.360
Pyrene	6.4	0.458

---

To understand how and where these parallel stacked molecules are situated around each other, we plot three dimensional spatial density functions. These represent an isosurface of the 3 dimensional distribution function of the centre of one molecule around another. The top panel of Fig. 4 shows the spatial density functions for all relative orientations of molecules in liquid phenanthrene and pyrene, plotted as the isosurface for the 15% most likely locations for a molecule in the first solvation shell. Compared to smaller aromatics, these solvation structures are more complex showing many diffuse lobes. A more detailed view of these diffuse lobes can be observed through xy and xz slices through the SDF's, as plotted in the SI Fig. S12–S15.

Fig. 4 Spatial density functions showing 15% most likely locations for all molecules in the first solvation shell in liquid phenanthrene (a) and pyrene (b), and for only parallel molecules (θ < 10°) for phenanthrene (c) and pyrene (d).

Following similar studies for benzene and naphthalene, we can enquire as to the location of parallel stacked molecules by plotting SDF's for molecules where the angle between the aromatic planes is less than 10°. This is shown in the bottom panel of Fig. 4. For phenanthrene the stacking is clearly heavily offset, with a preference for being above the concave “bay” section of the molecule and opposite this along the longer sides of the molecule. Conversely, for pyrene there is no clear offset in the stacking, with the preference being for face-to-face stacking. This is at odds with smaller aromatics benzene and naphthalene that do show a clear offset in parallel stacked molecules. In Fig. 5 we plot a 2d slice through the xz plane of the spatial density functions of parallel molecules, for all four aromatic liquids considered here (with x being the long axis in the plane of the molecule, and z being the axis normal to the aromatic plane). If we set aside phenanthrene (due to its decreased symmetry) there is a clear decrease in the offset of the stacking with increased aromatic size as one moves from benzene to naphthalene to pyrene (noting we get the same result here with the OPLS fixed charge used as the reference potentials, see SI Fig. S16). We can also check the relative orientation of the molecules by plotting the g(r,ϕ) where ϕ is the angle between vectors defining the long axis of the molecules (Fig. S11 in SI), there are clear peaks at 0 and 180° showing a preference of complete face-to-face alignment of the pyrene molecules.

Fig. 5 2D xz slices through the spatial density functions for benzene,¹² naphthalene,¹¹ phenanthrene and pyrene, where x is the long axis of the molecule and z is the vector normal to the aromatic plane.

In conclusion, we find that increasing the size and condensation of an aromatic from naphthalene to pyrene increases the amount of parallel stacking within the liquid structure. This can be simply rationalised by increased molecular size meaning that a larger number of atoms are pre-arranged for maximum dispersion and Coulomb attraction upon stacking. A similar result is also seen in the crystal structures, where benzene and naphthalene show herringbone structures of single molecules, pyrene shows a herringbone structure of stacked pairs offset along their long axis.¹⁷ Following from this, we find a clear trend that in the liquid state, stacking becomes distinctly overlapped with increasing aromatic size, with face-to-face aligned stacking observed as the highest probability stacking motif for pyrene. The reasons for this are less clear. Calculations of dimer structure show offset stacked or crossed structures as being energy minima.^7–9 However, the liquid state is more complex, with steric and/or packing considerations possibly favouring maximum overlap of the molecules. Finally the neutron data provide a robust experimental benchmark for comparing to simulation and theoretical predictions of intermolecular structure for these important and ubiquitous PAH systems.

Author contributions

TFH: conceptualisation, formal analysis, investigation, visualisation, writing – original draft, writing – reviewing and editing. NTS: investigation and writing – reviewing and editing. MPH: conceptualisation, investigation, writing – reviewing and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Collected data from NIMROD are available from “T. F. Headen, Polyaromatic Hydrocarbons: understanding liquid structure for benchmarking of atomistic and coarse-grained simulations, 2018, https://doi.org/10.5286/ISIS.E.RB1820596”, full files from Gudrun data reduction are available from https://github.com/tfheaden/neutrondata.

The data supporting this article have been included as part of the supplementary information (SI). The supplementary information contains a detailed description of the experimental and data refinement methods; fits to the data in real-space, with and without refinement and different potential models; full angular radial distribution functions; and selected 2D slices through the spatial density functions. See DOI: https://doi.org/10.1039/d5cc06373f.

Acknowledgements

The authors wish to thank Dr Camilla Di Mino for useful comments, and acknowledge the UK Engineering and Physical Sciences Research Council (grant EP/R013195/1) and the US National Science Foundation (grant CBET-1821762) for financial support, and UK Science and Technology Facilities Council (STFC) for award of beamtime (RB1820596).¹⁸

References

1. L. M. Salonen, M. Ellermann and F. Diederich, Angew. Chem., Int. Ed., 2011, 50, 4808–4842.
2. T. S. Totton, A. J. Misquitta and M. Kraft, Phys. Chem. Chem. Phys., 2012, 14, 4081.
3. A. G. Tielens, Annu. Rev. Astron. Astrophys., 2008, 46, 289–337.
4. M. Rapacioli, C. Joblin and P. Boissel, Astron. Astrophys., 2005, 429, 193–204.
5. O. C. Mullins, Annu. Rev. Anal. Chem., 2011, 4, 393–418.
6. T. M. Figueira-Duarte and K. Müllen, Chem. Rev., 2011, 111, 7260–7314.
7. L. Dontot, F. Spiegelman and M. Rapacioli, J. Phys. Chem. A, 2019, 123, 9531–9543.
8. E. M. Cabaleiro-Lago and J. Rodríguez-Otero, ACS Omega, 2018, 3, 9348–9359.
9. N. J. King and A. Brown, J. Phys. Chem. A, 2022, 126, 4931–4940.
10. H. Takeuchi, Comput. Theor. Chem., 2013, 1021, 84–90.
11. T. F. Headen, P. L. Cullen, R. Patel, A. Taylor and N. T. Skipper, Phys. Chem. Chem. Phys., 2018, 20, 2704–2715.
12. T. F. Headen, Mol. Phys., 2019, 117, 3329–3336.
13. S. E. Wheeler, J. Am. Chem. Soc., 2025, 147, 19738–19750.
14. K. Carter-Fenk and J. M. Herbert, Phys. Chem. Chem. Phys., 2020, 22, 24870–24886.
15. E. M. Cabaleiro-Lago, J. Rodríguez-Otero and S. A. Vázquez, Phys. Chem. Chem. Phys., 2022, 24, 8979–8991.
16. A. K. Soper, Phys. Rev. B, 2005, 72, 104204.
17. G. R. Desiraju, A. Gavezzotti and IUCr, Acta Crystallogr., Sect. B: Struct. Sci., 1989, 45, 473–482.
18. T. F. Headen, Polyaromatic Hydrocarbons: Understanding liquid structure for benchmarking of atomistic and coarse-grained simulations, 2018 DOI:10.5286/ISIS.E.RB1820596.

---