Combined Experimental and Computational Study of Polycyclic Aromatic Compound Aggregation: The Impact of Solvent Composition

Abstract

The aggregation of polycyclic aromatic compound (PAC) molecules is sensitive to the solvent they are dissolved or suspended in. By using both dynamic light scattering and diffusion-ordered nuclear magnetic resonance spectroscopy, in combination with molecular dynamics simulations, the effect of chemical structure on the aggregation of PACs in both aromatic and alkane solvents were systematically investigated. A suite of triphenylene-based PACs offers a robust platform to understand the driving forces of aggregation mechanism across both nanometer and micrometer scales. Both the configuration, either parallel or otherwise, and the arrangement, whether compact or loose, of molecules in their aggregates are determined by a fine balance of different interactions such as those between the polar groups, π–π interactions between the aromatic cores, steric hindrance induced by the side chains, and the degree of solvation. These results suggest that molecular architecture is the major factor in determining how the model compounds aggregate. The shift from aromatic to aliphatic solvent only slightly increases the likelihood of aggregation for the model compounds studied while subtle differences in molecular architecture can have a significant impact on the aggregation characteristics.

Keywords:

• Aggregation
• diffusion NMR
• DLS
• molecular dynamics

Introduction

The aggregation of polycyclic aromatic compounds (PACs) is determined by three factors: collision frequency between molecules, determined by their concentration, angle of interaction between molecules and the strength of molecular interactions, and balance of attractive and repulsive intermolecular forces. Attractive forces include van der Waals, π-π and polar interactions, hydrogen bonding, and metal coordination,^1–4 while repulsive forces are typically steric in origin, induced by any alkyl chains surrounding the aromatic core.⁵ The balance of these forces is highly sensitive to the solvent composition. When PACs are dispersed in an aromatic solvent, they are generally viewed as nanometer-sized aggregates^6–9: the strong affinity toward the solvent is likely attributed to the attraction between the aromatic solvent and the aromatic surface of the solute.¹⁰ However, introduction of n-alkane can alter the balance of solvent-solute interaction, facilitating the formation of micrometer-sized large aggregates.¹¹

This aggregation behavior, over two distinct length scales, has been observed in naturally occurring PACs such as asphaltenes. This broad family of compounds, one of the heaviest fractions of crude oil, are soluble in aromatic solvents, insoluble in aliphatic solvents, and known to aggregate readily upon changes in temperature, pressure, or solvent composition.^12–14 Aggregation, and subsequent precipitation, can occur during transportation or processing, resulting in blockage of pipelines and production facilities, fouling on equipment surface, and catalyst deactivation.^15–17 Studying the effects of solvent composition on asphaltene aggregation has significant importance on industrial practice such as blending.^18,19 To rationalize asphaltene aggregation in toluene solutions after n-heptane addition, a mechanism that is diffusion-limited below the critical micelle concentration but reaction-limited above has been proposed.^20–22 Other studies suggested that factors such as chemical structure, temperature, and the presence of co-solutes can influence the aggregation of asphaltenes.^23–29 For example, both the amount of asphaltene aggregation observed, and the rate at which it occurs, has been found to be inversely proportional to the chain length of the n-alkane added.³⁰

A number of experimental approaches, including dynamic light scattering (DLS),^31–33 nuclear magnetic resonance (NMR) diffusion-ordered spectroscopy (DOSY),^34,35 and small angle scattering, such as small angle x-ray scattering (SAXS)³⁶ and small angle neutron scattering (SANS),^37,38 have been used to investigate the aggregation process and to characterize the aggregates formed. For example, a combination of techniques was deployed by Espinat et al.³⁹ in studying the size of asphaltene aggregates in toluene as a function of temperature. They not only quantified the size of micelles in the region of tens of angstroms, but that different aggregation mechanisms could take place. In a separate work, Eyssautier et al.⁴⁰ used SANS and SAXS to investigate the structure of asphaltene nanoaggregates, describing the nanoaggregates as disks of total radius 32 Å, with 30% polydispersity, and a height of 6.7 Å. The reported sizes of asphaltene aggregates have ranged from a few nanometers to several microns, depending on the nature of the samples, the solvent quality, and other processing parameters.^40,41 Another approach that has been commonly used is molecular dynamics (MD) simulations where chemical structures of the same compounds were constructed and dispersed in different solvents to investigate the intermolecular interactions and molecular configurations.^42–47 The computational measurements offer a unique capability to quantitatively evaluate both aggregation distances and configurations between PAC molecules, and also those between solute and solvent molecules.

Due to the diverse range of chemical structures exhibited, the broad range of solvents that can be used, and the two distinctive length scales over which aggregation takes place, it is a challenging task to understand and control the solubility of natural asphaltenes, as work carried out previously has demonstrated. Instead, simple PAC molecules have been studied to unravel the fundamentals underpinning the aggregation mechanism of more complex PACs such as natural asphaltenes.^48,49 Either approach, using model compounds with controlled structure or investigating more disperse natural asphaltenes, has its limitations in revealing the aggregation mechanism of asphaltenes. In our recent work, a set of six model PAC molecules based on a triphenylene core were developed to study the aggregation mechanisms,⁵⁰ providing fundamental understanding of the same processes in other aromatic molecules including natural asphaltenes. Consistency in the core structures, and small, incremental changes in functionalized side chains, enables the effects of both functional groups and lengths of alkyl chains on aggregation to be evaluated separately.

In the work reported here, the effect of solvent composition, varying from purely aromatic (toluene) to purely aliphatic (n-heptane), on the aggregation characteristics of triphenylene-based PACs was investigated. DLS and diffusion-ordered nuclear magnetic resonance (DOSY NMR) measurements were used to monitor the hydrodynamic radii of any aggregates formed, while MD simulations were used to investigate the intermolecular interactions. This combined approach allows us to distinguish the role of the different solvent mixtures, differently sized alkyl chains, and different functional groups present and their effects on aggregation.

Materials and methods

Materials

Toluene and n-heptane of different grades (99.85%, Extra Dry-AcroSeal; 99+%, extra pure and deuterated) and PTFE membrane filters (100 nm pore size, Whatman) were purchased from Fisher Scientific (Loughborough, UK). Quartz optical cells (S High Precision Cell - Light Path 3 × 3 mm) were purchased from Hellman Analytics Q, and solvent resistant screw cap vials (7 and 15 mL) were purchased from Sigma Aldrich (Dorset, UK).

Triphenylene (TPN-C0) was purchased from Sigma Aldrich (Dorset, UK), from which two groups of model compounds were synthesized with either six alkoxy groups (TPN-C1, -C3, -C5, and -C10) or a single amide-linked chain (TPN-CN and -CNAcid) attached to the central triphenylene core. The chemical structures of TPN-C0 and all six model compounds are presented in Scheme 1. Details concerning the chemical synthesis can be found in a previous work.⁵⁰

Scheme 1. Chemical structures of the triphenylene (TPN) derivatives used in the present work.

Dynamic light scattering

A Zetasizer (Nano Series, Malvern) (λ = 632.8 nm), with a scattering angle of 173°, was used to carry out the DLS measurements, with data collected at fixed time intervals (immediately after n-heptane addition (0 h), 24 h, and 168 h). All samples were kept at constant temperature (22 ± 0.5 °C) and ambient pressure. Solvents were filtered three times with PTFE filters (100 nm pore size) prior to usage. The model compounds were stabilized in toluene for 24 h before n-heptane was introduced to form a mixture, the so-called heptol (4 parts toluene to 6 parts n-heptane), with a PAC concentration of 10 mg mL^‒1.

Analysis of the acquired autocorrelation functions for each sample was performed using the instrument software provided. Each datum is an averaged value of three independent samples, with each of them measured six times. Stokes–Einstein equation (Equation 1) was used to calculate the hydrodynamic radius of the aggregate, $r_{\mathrm{H}}$ (1) $$r_{\mathrm{H}}=\frac{kT}{6\pi \eta D}$$(1) where k is Boltzmann’s constant, T is the absolute temperature, and η is the viscosity of the solvent.^51,52

Diffusion-ordered NMR spectroscopy

Samples of triphenylene and the six different model PACs were prepared with a nominal concentration of 10 mg mL^‒1 in heptol-d (4/6 toluene-d₈:n-heptane-d₁₆ ratio), using tetramethylsilane (TMS) as reference (both Sigma Aldrich, Dorset, UK). Diffusion NMR data were collected at fixed time intervals (immediately after n-heptane addition (0 h), 24 h, and 168 h). All measurements were carried out on non-spinning on a 300 MHz Bruker Avance spectrometer, using a 5 mm PABBO BB-1H ZGRD probe equipped with a z gradient coil producing a maximum gradient of 36 G cm^‒1. All NMR measurements were performed at 295.15 K and used a double stimulated echo bipolar pulse pair sequence.⁵³ Diffusion coefficients and DOSY spectra were subsequently determined and processed using the DOSY Toolbox software package.⁵⁴ Further details on the diffusion NMR methodology used can be found in the Supporting Information.

The NMR experiments were used to confirm the presence of any nanoaggregates in the sample and estimate their size. For nanometer-sized particles, Equation (2)(2) $$D=\frac{kT}{6\pi \eta r_{\mathrm{H}}}\left(\frac{3\alpha }{2}+\frac{1}{1+\alpha }\right)$$(2) , (2) $$D=\frac{kT}{6\pi \eta r_{\mathrm{H}}}\left(\frac{3\alpha }{2}+\frac{1}{1+\alpha }\right)$$(2) the Gierer–Wirtz (G-W) modification of Equation (1),^55,56 where α is the ratio of solvent molecular weight to solute molar mass, provides more accurate predictions by taking the relative sizes of the solute and solvent into account. By combining the Gierer–Wirtz modification with appropriate assumptions about factors such as solute shape, flexibility, and solvation, a relatively simple method for predicting the diffusion coefficients of molecules based on their molecular weight can be obtained. From this predicted diffusion coefficient, it is possible to estimate a hydrodynamic diameter for the dis-aggregated molecules, to be compared with the experimentally acquired data.

Molecular dynamics simulations

Molecular dynamics (MD) simulations were carried out using the GROMACS 4.6.5 software package.^5,42,57,58 Starting from an initial energy minimization procedure using the steepest descent algorithm (10⁴ iterations, convergence criterion of 25 kJ/mol/nm and initial step size of 0.01 nm),⁵⁹ short sequential simulations (0.1 and 3.0 ns, respectively) in the NVT and NPT ensembles were used to relax the systems. A Berendsen barostat algorithm and a stochastic velocity rescale⁶⁰ were adopted to maintain 1 atm pressure and 298 K respectively. At each stage, equilibration was confirmed using the algorithm developed by Chodera⁶¹ and implemented in the pymbar python library⁶² (data not shown here). Once equilibrated, production simulations lasting 20 ns in the case of one TPN molecule and 100 ns in the case of multiple were performed with an integration time of 2 fs using a Parrinello-Rahman barostat and the stochastic velocity rescale thermostat. Throughout these simulations, the dimensions of the simulation box were 5 nm (±0.3 nm) in each direction.

Bonding interactions were modeled using the OPLS/AA forcefield, as reported in previous studies.^44,63 Molecular structures of the six different compounds and the triphenylene precursor were constructed using the protocols established in a previous study of this group.⁶⁴ Non-bonding Van der Waals interactions were modeled using the 12-6 form of the Lennard–Jones functional, and partial atomic charges were modeled using the particle-mesh Ewald algorithm. Parameters for Van der Waals and electrostatic nonbonding interactions were again taken from the OPLS-AA forcefield, with partial charges modified to ensure charge-balancing. A cutoff of 1.0 nm for nonbonding interactions was used with potential is shifted such that the value was zero at the cutoff, and three dimensional periodic boundary conditions were used for all equilibration and production simulations.⁶⁵ The effectiveness of such fixed-charge forcefields for accurately modeling aggregation of extended aromatic systems is subject to debate.^66,67 For the purposes of this study, these interactions should provide insight into the trends observed in the experimental system even if quantitative details (e.g. aggregate size, π–π interaction distances) are not necessarily reproduced.

Solute-solvent affinity was studied by placing a single PAC molecule in a simulation box filled with 700 solvent molecules (350 n-heptane and 350 toluene), from which a python package MDAnalysis^68,69 was used to evaluate the solvation state around the molecule. Shells were set at distances of 0.5, 1.0, and 1.5 nm from the outermost atom in the molecular structure of the PAC. The number of toluene and n-heptane molecules in each shell was recorded every 1 ns over a total 100 ns of the simulation time, from which the fraction of toluene molecules was recorded.

To investigate the intermolecular interactions between solutes, seven identical molecules were placed into (1) toluene; (2) n-heptane simulation boxes (each containing 700 solvent molecules); and (3) a heptol simulation box (containing 350 toluene/350 n-heptane molecules) at equidistant positions to examine the intermolecular forces. Potential of mean force (PMF) analysis was used to determine the averaged distance between the centers of mass (COMs) of molecules during the simulation, based on the data extracted from the radial distribution function of GROMACS. The positions and intensities of the distribution peaks acquired offer provide the averaged distance between the molecules. The data obtained has been normalized against the interaction recorded at 2.4 nm. In the case of the heptol samples, additional MD analysis provided a more detailed insight into the molecular configurations adopted during aggregation. The nature of aggregation, i.e. π-stacking, ordered, dis-ordered, was determined by capturing both the intermolecular distances and the angle between aromatic cores during the simulations. This set of data was only extracted for the heptol boxes and is presented as histograms to support the PMF potential energy curves.

It is worth noting that a 50/50 molar ratio was used in the MD simulations to match the results between the solvent affinity (1 model compound molecule in 700 solvent molecules) and solute-solute (7 model compound molecules in 700 solvent molecules). However, 60% of the toluene was replaced with n-heptane in the DLS and NMR experiments based on the consideration that most studies reported that n-alkane solvent effects are exhibited when more than 50% of the solvent volume is n-heptane.^6,70,71

Results

Dynamic light scattering of PACs

The hydrodynamic diameters of PAC aggregates in heptol solution were measured as a function of time, using a protocol established previously, where the aggregation of the same set of model compounds was studied in toluene.⁵⁰ The averaged hydrodynamic diameters acquired over a course of 168 h (7 days) are presented in Figure 1.

Figure 1. Averaged hydrodynamic diameters of seven model PACs dispersed in heptol (40% toluene, 60% n-heptane), at 10 mg mL^‒1. Measurements of diameters smaller than 1 nm, typically beneath the detection limit of the instrument, are marked with an asterisk.

It was found that TPN-C0, the unaltered triphenylene, presented a hydrodynamic diameter below 1 nm, beneath the typical detection limit of the instrument, in heptol solution with no noticeable change over 168 h, consistent with the previous finding that TPN-C0 does not aggregate in toluene⁵⁰ and confirming that adding n-heptane does not prompt aggregation of the TPN-C0. In contrast, the compounds with alkyl branches, including TPN-C1, -C3, -C5, and -C10, were all found to form large aggregates, hundreds of nanometers in size, upon the addition of n-heptane. TPN-C3 and TPN-C5 were subsequently found to dissociate into smaller aggregates over an extended period, different aggregation kinetics to those observed in toluene. Compared with the observation made in pure toluene,⁵⁰ whereby the aggregates of TPN-C3 and -C5 dissociated after 24 h, it took much longer for these aggregates to dissociate in heptol. These results highlight the influence of solvent quality on both aggregation and any later dissociation.

The PACs carrying polar functional groups, TPN-CN and TPN-CNAcid, were found to form large aggregates in heptol with hydrodynamic radii of ca. 600 and 1400 nm, respectively. After 168 h, TPN-CNAcid exhibited a significant reduction in particle size, to ca. 300 nm, while TPN-CN aggregates expanded to over 1000 nm. Such variations in particle size, evidenced by the large standard errors, imply the likelihood of a complex aggregation mechanism, resulting in a nonuniform size distributions, even after 168 h. As TPN-CNAcid is considered of greater polarity than TPN-CN, the differences observed at 168 h could be attributed to the precipitation of the clusters, consistent with previous studies that suggest that natural asphaltenes of different origins exhibit different precipitation characteristics.^72–74

Diffusion NMR of PACs

Diffusion NMR data of the set of model compounds were acquired as a function of time to supplement the light scattering measurements. The use of a mixed solvent, heptol, complicated the analysis of the diffusion NMR data because additional signals are introduced to the ¹H spectra. Figure 2 shows a representative ¹H DOSY spectrum of TPN-C5 in heptol-d. The ¹H NMR spectrum of TPN-C5 shows the distinct chemical regions of the model compounds: aromatic groups at ca. 8 ppm, methoxy-groups at ca. 4 ppm, and alkyl chains at <2 ppm. However, the presence of n-heptane introduces three additional peaks in the alkyl region, resulting in overlapped signals in the ¹H spectra and subsequent compromise diffusion coefficients in the DOSY spectra. Aromatic and methoxy peaks remain well-resolved and can be used to reliably estimate the diffusion coefficients for TPN-C0, TPN-C3, TPN-C5, and TPN-C10. However, it was not possible to obtain diffusion data of satisfactory quality for all the compounds studied. TPN-C1, TPN-CN, and TPN-CNAcid were found to form very large aggregates in both toluene and heptol solutions. The ¹H NMR spectra of these samples were dominated by the solvent signals, while the peaks attributed to the PAC molecules themselves were very low in intensity, with broad signals in both the aromatic and aliphatic regions, resulting in unreliable estimation of diffusion coefficients. No diffusion data from these three species was used in this analysis. All ¹H and DOSY spectra used in this study can be found in the Supporting Information.

Figure 2. ¹H DOSY spectrum of 10 mg mL⁻¹ solution of TPN-C5 in heptol-d. Residual toluene and heptane signals are indicated by T and H respectively.

Figure 3 compares the hydrodynamic diameters estimated from the experimentally acquired diffusion coefficients for TPN-C0, TPN-C3, TPN-C5, and TPN-C10 with those predicted from the molecular weights of model compounds, using Equation (2)(2) $$D=\frac{kT}{6\pi \eta r_{\mathrm{H}}}\left(\frac{3\alpha }{2}+\frac{1}{1+\alpha }\right)$$(2) . As the SEGWE prediction yields the expected hydrodynamic diameter for a single molecule, any observed differences in diffusivity and, implicitly, hydrodynamic size can be attributed to the existence of nanoaggregates. Since TPN-C0 exhibited diffusion coefficients very similar to those predicted by Equation (2)(2) $$D=\frac{kT}{6\pi \eta r_{\mathrm{H}}}\left(\frac{3\alpha }{2}+\frac{1}{1+\alpha }\right)$$(2) , it does not form stable nanoaggregates, consistent with the previous DLS results. The DOSY-NMR results for TPN-C3, TPN-C5, and TPN-C10, however, each indicate that the species in solution move slower than expected and, therefore, confirm the presence of nanoaggregates. The observed diffusion coefficients in a DOSY spectrum are a weighted average of the diffusion coefficients of all the species contributing to a given signal in the NMR spectrum. These data, therefore, also indicate that any aggregates formed are likely in dynamic equilibrium with single species. However, these aggregates are markedly different to, and much smaller than, the micrometer sized clusters indicated by the DLS results in Figure 1. The solvent diffusion coefficients remain the same for all samples, indicating that any aggregates formed do not contain solvent molecules within the aggregates. Any discrepancy between results generated by DOSY NMR and DLS is not entirely surprising because the two techniques capture different fractions of the aggregates distribution: nanoaggregates are best detected by NMR, while large clusters can only be observed by DLS. It also highlights the importance of covering multiple length scales when studying the aggregation mechanisms of PACs.

Figure 3. Experimentally-acquired hydrodynamic diameters for TPN-C0, -C3, -C5, and -C10 in heptol-d solution at 0 (black), 24 (light grey), and 168 h (dark grey) after sample preparation. Predicted hydrodynamic diameters for single model compound molecules in heptol-d indicated by horizontal black bar.

Molecular dynamics simulations: Solute–solvent interactions

MD simulations were deployed to quantify the interactions between model PACs and both solvents present in the heptol solution. This data can be presented as a series of snapshots of the simulations (Figure 4) that offer a visual guide to the relative amounts of toluene (grey) and heptane (green) surrounding a PAC molecule, in this case, TPN-C1. Toluene molecules appear to dominate in regions close to the molecule. Further away, the difference in affinity reduces and equal numbers of the two solvents are present.

Figure 4. Set of snapshots of TPN-C1 in 50/50 heptol mixture, produced using VMD 1.9.3, with toluene molecules (gray) and heptane molecules (green) shown within (a) 0.5, (b) 1, and (c) 1.5 nm of the TPN molecule.

Such data can also be presented as histograms (Figure 5) whereby the percentage of toluene or n-heptane within each shell indicates the relative affinity of solvent molecules toward the model compound, with 50% indicating an equal affinity between solvents. This analysis methodology was used previously⁷⁵ to evaluate molecular interactions with aliphatic and aromatic hydrocarbons.

Figure 5. Percentage of toluene within successive shells 0.5, 1.0, and 1.5 nm from the aromatic surface of PAC molecules in a heptol simulation box over 100 ns simulation time. Horizontal dashed line indicates an equal toluene/n-heptane distribution. Error bars are 95% confidence interval around the mean (n = 25 samples).

The unfunctionalized model compound (TPN-C0) showed a preference for n-heptane at all three distances sampled. However, variation in chemical structure, such as functionalizing the fused aromatic rings with even only methoxy groups, influences the preference toward toluene. For example, TPN-C1 showed the strongest affinity of all the model compounds toward toluene, with approximately 70% toluene within the 0.5 nm shell, decaying to 50% after 1.0 nm. The unexpected behavior of simulated TPN-C0 is likely caused by the lack of alkoxy substitution and consequent electron-withdrawing effect on the aromatic core. Alkoxy substitution leads to a relative electron deficiency on the TPN core, represented in the forcefield as partial positive atomic charges, and increased the affinity toward electron-rich toluene.

As the aliphatic chain length was increased in TPN-C3, -C5, and -C10, the strong preference for toluene was reduced. It is probable that an increased chain length of aliphatic branches enables closer association between n-heptane and the PAC, while at the same time repels the toluene. Longer aliphatic chains also limit the solute-toluene interaction to only the aromatic area of the model compound, contributing to the reduced presence of toluene close to the aromatic core.

Neither TPN-CN nor TPN-CNAcid (Figure 5) showed any noticeable preference between the aromatic and alkane solvent, both of which are less polar than the amide and acid groups present. To further evaluate the effect of the polar functional groups on solvent preference, the analysis was repeated for the surface of the aromatic core and the aliphatic tails of both TPN-CN and TPN-CNAcid (Figure 6). This indicates that the tail of TPN-CNAcid has a stronger affinity toward toluene than that of TPN-CN.

Figure 6. The percentage of toluene within successive shells 0.5, 1.0, and 1.5 nm from (left) the aromatic surface and (right) the aliphatic tail group of PAC molecules TPN-C0, -CN, and -CNAcid, suspended in a heptol simulation box over 100 ns simulation time. Horizontal dashed line indicates an equal toluene/n-heptane distribution. Error bars are 95% confidence interval around the mean (n = 25 samples).

Molecular dynamics simulations: Intermolecular interactions between PACs

To evaluate the effect of solvent on the intermolecular interactions between the PACs, MD simulations were carried out in three solvent boxes. PMF analysis was performed for all seven compounds, estimating the magnitude and nature of the intermolecular interactions between the PACs as a function of the solvents used. The potential of mean force curves of TPN-C3, C5, and C10 all feature a series of evenly spaced potential energy minima. These reveal the presence of π-stacking between individual model compound molecules and confirm the tendency to form nanometer-sized aggregates. Figure 7 shows a representative PMF analysis extracted from the simulation data of TPN-C3 in pure toluene (dark brown), 50/50 toluene/n-heptane (rust brown), and pure n-heptane (red) respectively. In addition to the distance and the magnitude of the interactions (Figure 7a), statistical analysis of the distance between molecules (Figure 7b) and the angles between the aromatic cores of PAC molecules (Figure 7c) is presented. Figure 7b shows explicit dimerization in the mixed solvent system, in accordance with the PMF, which is supported by the angle histogram (Figure 7c), which shows the TPN-C3 molecules aggregated in a parallel configuration, i.e. through π-stacking.

Figure 7. (a) Potential of mean forces between the centers-of-mass of TPN-C3 in 700 toluene molecules (dark brown), heptol (350 toluene molecules and 350 n-heptane) (rust brown) and 700 n-heptane molecules (red) over 100 ns simulation time. Data were normalized to 0 at 2.4 nm. (b) Histogram of aggregation distance between the centers of mass of TPN-C3 molecules in heptol during the simulation. (c) Histogram of angles between aromatic cores of TPN-C3 molecules in heptol during the simulation. Angle data was only recorded when aromatic groups are within 1.5 nm of each other.

Different PMFs were observed for the two PACs with polar groups, TPN-CN and TPN-CNAcid. Both compounds exhibited broad PMFs, with a peak at ca. 0.5 nm and significant intensity at 1.0 nm and beyond. As the amount of n-heptane in solution is increased, while the major peaks did not change position, the depth of the interaction increased and shoulder peaks developed at distances >0.5 nm. A representative dataset for TPN-CNAcid is shown in Figure 8.

Figure 8. (a) Potential of mean forces between the centers-of-mass of TPN-CNAcid in 700 toluene molecules (dark brown), heptol (350 toluene molecules and 350 n-heptane) (rust brown) and 700 n-heptane molecules (red) over 100 ns simulation time. Data was normalized to 0 at 2.4 nm. (b) Histogram of aggregation distance between the centers of mass of TPN-CNAcid molecules in heptol during the simulation. (c) Histogram of angles between aromatic cores of TPN- CNAcid molecules in heptol during the simulation. Angle data was only recorded when aromatic groups are within 1.5 nm of each other.

This analysis was also performed for TPN-C0, -C1, -C5, -C10, and -CN, for which a comprehensive set of figures based on MD simulations, including both potential of mean force analysis and radial distribution functions, can be found in the Supporting Information. The key features of this dataset, i.e. distances between molecules and the magnitudes of potential energy minima, are summarized in Figures 9 and 10.

Figure 9. Histograms of the interaction strength between PAC molecules in simulation boxes of (a) toluene, (b) heptol, and (c) n-heptane, normalized against that recorded at 2.4 nm. More than one interaction minima were observed for some compounds, indicating the existence of either multiple aggregation configurations or multiple molecules in one stack.

Figure 10. Histograms of the molecular distances corresponding to the potential energy minima between PAC molecules in simulation boxes of (a) toluene, (b) heptol, and (c) n-heptane.

The magnitude of the background thermal energy (1 kT) of all molecules examined was approximately 2.5 kJ mol^‒1 at 298 K. Any observed interaction lower than 2.5 kJ mol^‒1 will be insignificant. The interaction strengths between TPN-C0 molecules are smaller than this benchmark for all solvents used, toluene, heptol, and heptane, indicating that TPN-C0 molecules show little attraction to each other in any of the solvents used. This is consistent with experimental NMR and DLS data that indicates there was barely any aggregation observed for this molecular structure. The histogram (Figure SI.2.1) shows that the distance between TPN-C0 molecules followed a unimodal distribution centered around 3.0 nm. In the few cases when TPN-C0 molecules came close to one another, no preferred orientation, e.g. parallel or perpendicular stacking, was observed.

In contrast, TPN-C1 molecules in toluene exhibited strong interactions over a broad range of intermolecular distances between 0.7 and 2.0 nm, with a minimum interaction energy of −12.0 kJ mol^‒1 observed at 0.85 nm (Figure SI.2.2). This leads to a significant aggregation observed in the simulation carried out, in agreement with previous MD studies of PAC aggregation^42–44 and echoes the results of the experimental measurements (Figures 1 and 3). Both the potential of mean force curve and the corresponding statistical analysis suggest that the parallel configuration was not favored by TPN-C1 molecules in toluene. However, once toluene was replaced by n-heptane, the interaction energies between the TPN-C1 molecules became stronger while the fine structure reveals several narrower minima. This is highly likely due to the changed configuration of TPN-C1 molecules in the aggregates, from loose to compact, as they become surrounded by an excessive amount of n-heptane. This is consistent with both the solvent affinity simulation that implies an increased affinity toward toluene and the light scattering data that shows a reduced particle size as n-heptane is added.

Series of narrow energy minima, indicative of parallel configurations, were observed from the PMF analysis of TPN-C3 (−11.5 kJ mol^‒1 at 0.39 nm (Figure SI.2.3)), TPN-C5 (−12.2 kJ mol^‒1 at 0.42 nm (Figure SI.2.4)), and TPN-C10 (−9.0 kJ mol^‒1 at 0.45 nm (Figure SI.2.5)), as presented in Figure 9. The first energy minimum for each of these molecules presents a shoulder, which is likely caused by the steric effects induced by the alkoxy tails presented around the aromatic core. The PMF curves also show secondary and tertiary interaction minima, indicating the parallel stacking of multiple molecules. The histograms of molecular orientation/angles for TPN-C3, -C5, and -C10 present a strong preference toward parallel configurations of the aromatic centers, supporting the hypothesis that the PACs form π-stacked aggregates. The magnitude of the energy minima increased both as the length of the alkyl branch chains was increased from propyl to pentyl and when n-heptane was added. For TPN-C10 in heptane, however, the fine characteristics of successive energy minima were no longer observed, replaced with a single broad interaction peak, which suggests that the favorable configurations were much less well-defined, and the resulting aggregates or clusters were not structurally ordered.

Both TPN-CN and TPN-CNAcid exhibited broad, rather featureless, potential energy curves (Figures SI.2.6 and SI.2.7), which suggests that they are prone to forming aggregates with a less uniform structure than those of the alkoxy-modified compounds. A major potential energy minimum was found between 0.3 and 0.4 nm, the characteristic feature for π-stacking. However, the much broader interaction peak observed, supported by the histogram data on angles of interaction, indicates a much greater range of possible configurations, including both head-to-tail and head-to-head arrangements, that TPN-CN and TPN-Acid molecules adopt when forming aggregates. The addition of n-heptane increased both the width and magnitude of the energy minima, suggesting a higher tendency for these model compounds to aggregate in alkane solutions.

Discussion

Effect of chain length on aggregation in mixed solvents

TPN-C0: As the benchmark compound, experimental hydrodynamic diameters determined by both DLS (Figure 1) and diffusion NMR (Figure 3) indicate that the presence of 60% n-heptane was not sufficient to facilitate aggregation. MD simulations of TPN-C0 confirm that aggregation occurs neither in a mixed heptol solvent nor in pure n-heptane. PMF data (Figures 9 and SI.2.1) does show a minor attractive interaction between TPN-C0 molecules, albeit lower than the background thermal energy of the systems studied. The finding is supported by the angle and distance histograms captured, showing a low occurrence of collisions and an absence of any preferred angle of interaction. Previous studies proposed that whether π-stacking is or is not the major driving mechanism for aggregation in natural asphaltenes depends on the size/number of the aromatic rings. With only a limited number of aromatic rings, TPN-C0 likely does not have an aromatic surface that is large enough to facilitate π-stacking.^4,69,76–80

Tpn-C1: DLS data suggest that the addition of n-heptane results in aggregates with slightly smaller hydrodynamic diameters than those found in pure toluene.⁵⁰ Such an aggregation mechanism is consistent with previous studies where the addition of the alkoxy groups surrounding the triphenylene core was found to facilitate the aggregation between the aromatic molecules.^81,82 Solvent affinity simulations (Figure 5) suggest that TPN-C1 molecules experienced a far greater attraction toward toluene than n-heptane. Interactions between the methoxy groups appear to counteract the aromatic-aromatic interactions between the relatively small aromatic surface of TPN-C1. No angle dependence of TPN-C1 aggregation was observed (Figure SI.2.2), indicating that there is very little π-stacking between the TPN-C1 molecules and that random configurations are preferred. This is likely because the methoxy groups present in TPN-C1 are not long enough to impose a stacking configuration that would not only limit the size of the aggregates but also stabilize the self-associated assembles.

TPN-C3 and -C5: Once the length of the alkoxy branch was further extended to propyl and pentyl, the PAC molecules showed less sensitivity between aromatic and n-alkane solvents than TPN-C1 (Figure 5), and the change in solvent did not facilitate aggregation. The hydrodynamic diameters of TPN-C3 and -C5 recorded by DLS (Figure 1) and DOSY (Figure 3) in heptol show the existence of both large clusters and nanoaggregates. Multiple minima, spaced by a consistent intermolecular distance of ca. 0.4 nm, were observed in the PMF data of TPN-C3 and -C5 in both heptol and pure n-heptane (Figures 7, SI.2.3, and SI.2.4), evidence for an explicit preference toward π-stacking⁷⁹ in agreement with previous experimental results confirming that TPN-C3 forms columnar crystalline phases while TPN-C1 cannot due to insufficient steric ordering.⁸³ Both the close range and low intermolecular angle recorded from the PMF simulations of TPN-C3 and -C5 acquired in heptol confirm the π-stacking configuration. The addition of heptane increases the magnitude of the potential energy minimum.

TPN-c10: Despite one of the six propyl branches being replaced by a decyl chain, TPN-C10 exhibited very similar characteristics to that of both TPN-C3 and TPN-C5. DLS results show that TPN-C10 molecules formed large clusters in heptol which remained stable for an extended period of time. NMR experiments captured the coexistence of nanoaggregates, likely consisting of only a few molecules stacked closely together. All of the MD simulations of TPN-C3, TPN-C5, and TPN-C10 (Figures SI.2.3–SI.2.5) in toluene and heptol show similar behavior: a series of potential energy minima of decreasing depth spaced by ca. 0.4 nm, a characteristic of π-stacking. Therefore, the increased hydrodynamic diameter of TPN-C10 aggregates that remained stable for 168 h can be attributed to the single increased chain length, the only difference between TPN-C3 and TPN-C10. The angle and distance histogram results of TPN-C10 suggest that planar configurations within the aggregates remain the preferred option, which is consistent with previous studies that have shown that longer tails limit potential configurations.⁵ The finding established from TPN-C10 indicates that non-centrosymmetric molecular structures with long side chains could have an increased affinity toward n-alkane solvent species, as would be expected of a model compound that replicates asphaltenes that are soluble in toluene but insoluble in n-heptane.^21,84,85

Solvophobic effects were only observed for TPN-C0 and TPC-C1. For TPN-C3, TPN-C5, and TPN-C10, the formation of nanoaggregates occurs in a similar manner in n-heptane as it does in toluene, with the surrounding alkyl chains imposing a π-stacking configuration. However, a slow dissociation of the large clusters initially formed was recorded by DLS measurements. It is very probable that once the alkyl tail-groups are long enough to favor a parallel (π-stacked) configuration, the nature of the solvent has a limited impact on the aggregation mechanism.^81–83 Changes in solvent composition appear to change how aggregation progresses, e.g. forming clusters with a great aggregation number, rather than initiating the aggregation.

Effect of amide and acid groups on aggregation in mixed solvent

TPN-CN and -CNAcid: The presence of polar functional groups, such as amides (TPN-CN) and acids (TPN-CNAcid), significantly increases the likelihood of aggregation of aromatic hydrocarbons in toluene.⁵⁰ While the aggregates formed in heptol were too large to measure using diffusion NMR, DLS results (Figure 1) show the presence of larger aggregates than in pure toluene for both TPN-CN and TPN-CNAcid, with a notable difference between their kinetics: the sizes of TPN-CN aggregates appeared to progressively increase, while those of TPN-CNAcid aggregates decreased over the same period.

Solvent affinity simulations (Figure 5) suggest that neither TPN-CN nor TPN-CNAcid exhibits an excessively strong attraction or repulsion toward either toluene or n-heptane. However, when studying the solvent affinity of the aliphatic moieties only (Figure 6), the acid and amide groups exhibit a stronger affinity toward toluene, indicating that the greater electron density of the acid and amide groups attract the aromatic solvent. There was no evidence of π-stacking in the PMF results for either TPN-CN or TPN-CNAcid (Figures 9, 10, SI.2.6, and SI.2.7) in any of the solvent mixtures. The histograms obtained from the MD simulations in heptol indicate no preferred angles of interactions and the configurations observed in MD simulations showed both tail-to-head and head-to-head formations. These results suggest that the primary driving force for aggregation of TPN-CN and TPN-CNAcid lies in the functional groups, rather than the aromatic cores.

The PMF analyses of TPN-CN and TPN-CNAcid display similar changes as the n-heptane was introduced to toluene. In pure n-heptane, there is a significant increase in both magnitude and width of the potential energy minimum, suggesting that the added alkane molecules facilitated greater interaction between the aromatic compounds, which is consistent with what has been observed in natural asphaltenes.⁶ The potential energy curve for TPN-CN data exhibits only a small change in the heptol solvent but a larger one when toluene is completely replaced by n-heptane. The aggregation of TPN-CN appears to be disordered in all solvents studied. While the PMF analysis does show a deep interaction centered at 0.5 nm, this is very broad and has no significant peaks in the angle histogram, consistent with a previous study.⁷⁹ By contrast, the PMFs of TPN-CNAcid shows a gradual change as n-heptane replaces toluene. The interaction peak in the PMF remains very broad, implying a number of possible configurations, nevertheless the TPN–TPN interaction energy, indicating its propensity to aggregate, increases.

Direct comparison between TPN-CN and TPN-CNAcid supports the hypothesis that acid groups play a major role in the aggregation of PACs. They not only facilitate aggregation, but also lead to significantly different configurations, which subsequently alter the macroscopic properties and behavior of the aggregates formed. This is consistent with previous studies whereby model PACs containing acid functional groups were observed to aggregate at oil/water interfaces; the aromatic cores subsequently driving a reconfiguration process that aligned the molecules in parallel configurations.^48,86,87 In these studies, the increased polar attractions proved the substantial drive to aggregation rather than the electron clouds present on the aromatic surface.^88,89 However, once aggregated, the π-stacking interactions facilitate the re-configuration of the aggregates into the most energetically favorable configuration, parallel, due to the aromatic cores of the molecules.⁹⁰ This conclusion is consistent with recently reported experimental observations concerning asphaltenes that are either strongly or weakly interfacially active.^91,92

Conclusions

A combination of experimental and computational approaches was deployed to systematically investigate the effect of solvent composition on the aggregation characteristics and kinetics of a set of model PACs. Light scattering and diffusion NMR data indicated that the hydrodynamic diameters of the aggregates cover a broad range, from a few nm to several μm, of which the molecular details were studied by MD simulations.

The aggregation of polycyclic aromatic hydrocarbons is determined by a balance between several different interactions: (1) π–π interactions between the aromatic cores, (2) interactions between any polar groups present, (3) steric hindrance caused by the side chains, and (4) the degree of solvation by toluene. The results acquired in the present work suggest that the structures of the compounds are the major factors in how the model compounds aggregate. The presence of the polar groups, both amide and acid, facilitates the formation of macroaggregates. While polar attractions overshadow the molecular attractions driven by van der Waals interactions, they are generally less sensitive to the solvent used. Side chains facilitate the formation of π-stacked nanoaggregates, however, this behavior was only observed when the alkyl groups are suitably long. It appears to require the presence of a long, possibly asymmetric, alkyl chain for the addition of n-heptane to induce the formation of macroaggregates in this sub-set of the model compounds. For the compounds that initially form large aggregates in toluene, the change in solvent promotes the development of larger aggregates still.

Naturally occurring PACs, such as asphaltenes, possess a wide range of aromatic cores, side chains of different lengths and polar groups, and are much more complex than the molecules investigated here. However, the methodology and the underpinning knowledge established through this work could be readily applied to polycyclic aromatic hydrocarbons of increasing complexity, with variations in the nature of the polar groups present, any alkyl chains attached and their dispersity, or the size of the aromatic cores.

Author contributions

The manuscript was written through contributions of all authors and all authors have given approval to the final version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by the Engineering and Physical Science Research Council under Grant EP/P007864/1. The diffusion NMR experiments were funded by the RSC Research Fund under Grant RF17-3528.