Aromatic sulfonation with sulfur trioxide: mechanism and kinetic model

Abstract

Electrophilic aromatic sulfonation of benzene with sulfur trioxide is studied with ab initio molecular dynamics simulations in gas phase, and in explicit noncomplexing (CCl₃F) and complexing (CH₃NO₂) solvent models. We investigate different possible reaction pathways, the number of SO₃ molecules participating in the reaction, and the influence of the solvent. Our simulations confirm the existence of a low-energy concerted pathway with formation of a cyclic transition state with two SO₃ molecules. Based on the simulation results, we propose a sequence of elementary reaction steps and a kinetic model compatible with experimental data. Furthermore, a new alternative reaction pathway is proposed in complexing solvent, involving two SO₃ and one CH₃NO₂.

---

Introduction

Aromatic sulfonation is a very important chemical transformation of organic compounds.^1,2 It belongs to the well-known class of electrophilic aromatic substitution (S_EAr) reactions alongside nitration, halogenation, acylation and alkylation. Sulfonation is a key reaction step in large industrial applications including pharmaceuticals, detergents, surfactants, dyes, and pesticides.^3–5 The most commonly used sulfonating agents are sulfur trioxide (SO₃), oleum, sulfuric acid and chlorosulfuric acid.⁶

Traditionally, S_EAr reactions are explained by a two-step S_E2 or arenium ion mechanism.⁷ In the first and rate-determining step, an electrophile attacks the electron-rich aromatic ring to form a σ-complex or wheland or arenium ion intermediate, which is stabilized by mesomery. Aromaticity is restored in the second step by elimination of H⁺. Recently, alternative pathways to the arenium mechanism have been proposed by Schleyer and coworkers in various electrophilic aromatic substitution reactions⁸ including halogenation,^9,10 nitration,¹¹ and sulfonation.¹² These studies highlight the diversity in reaction mechanisms for S_EAr and their dependence on substrate and reaction conditions.

Based on kinetic experiments on chlorobenzene and 1,4-dichlorobenzene, Cerfontain and coworkers proposed a three-step kinetic scheme for the sulfonation of arenes with SO₃ in aprotic solvents (Scheme 1).^13–16 The first step is a reversible reaction wherein a σ-complex is formed between the arene and a SO₃ molecule. In the second step a second SO₃ reversibly binds to the first SO₃. Finally, proton transfer from the arene to the second SO₃ in the third step restores aromaticity and drives the reaction toward arenepyrosulfonic acid (ArS₂O₆H), which can be readily hydrolyzed in aqueous media. In apolar noncomplexing CCl₃F solvent sulfonation kinetics are first order in SO₃, and thus step 1 was thought to be rate-limiting.¹⁷ In contrast, in polar SO₃-complexing¹⁸ CH₃NO₂ the rate is second order in SO₃, and step 2 was taken as the rate-limiting step.¹⁷

Scheme 1 Kinetic model proposed by Cerfontain and coworkers for sulfonation with SO₃.

If less than two equivalents of SO₃ are used per mole arene, two sulfonation stages can be distinguished. A fast primary stage purportedly proceeds via the reaction steps 1–3 (Scheme 1). A much slower secondary stage was proposed to proceed through a reaction of arenepyrosulfonic acid with arene, forming two molecules of arenesulfonic acid (ArSO₃H, step 4).^13,19 This study, however, is focused on the primary sulfonation stage.

The mechanism proposed by Cerfontain and coworkers involving a σ-complex intermediate with one SO₃ has been challenged by recent theoretical studies. Morley et al.^20,21 performed quantum chemical calculations at the Hartree–Fock level on the sulfonation of toluene with SO₃ in gas phase, and concluded that formation of a toluene–SO₃ σ-complex was unlikely due to the high energy change required. Instead, they proposed the initial formation of a toluene⋯SO₃⋯SO₃ π-complex followed by a toluene–S₂O₆ σ-complex with almost the same energy, and a cyclic proton rearrangement yielding toluenepyrosulfonic acid.

Schleyer and coworkers studied the sulfonation of several arenes with static density functional and SCS-MP2 calculations in implicit solvent.¹² The sulfonation of benzene, 1,4-dichlorobenzene, toluene, and naphthalene, in gas phase and in apolar solvent were found to proceed via a concerted pathway involving two SO₃ molecules forming a cyclic σ-complex transition state (TS), without intermediate (Scheme 2). In CH₃NO₂, the cyclic σ-complex became an intermediate state, but with low stability. With only one SO₃, no intermediate σ-complex was formed in both solvents, and very high energy barriers were needed for sulfonation.

Scheme 2 Concerted sulfonation mechanisms in gas phase or apolar solvent with 2 SO₃ molecules according to Schleyer and coworkers,¹² and with 1 SO₃ + H₂SO₄ catalyst according to Morkovnik and Akopova.²³

The absence of an intramolecular primary hydrogen kinetic isotope effect (KIE)¹⁷ led Cerfontain and coworkers to conclude that the reaction step involving proton transfer (step 3 in Scheme 1) is not rate-limiting. However, the KIE may also be small if the C–H bond is only partially broken at the TS.²² In the cyclic transition state structures involving two SO₃ molecules, as calculated by Schleyer and coworkers,¹² the C–H bond is only slightly elongated at the TS in both CCl₃F (0.08 Å) and CH₃NO₂ (0.10 Å at the second TS). The small bond elongations are consistent with ratios k_H/k_D = 1.2–1.3, suggesting that step 3 may indeed be rate-limiting.

Importantly however, it is not clear how the participation of two SO₃ molecules, as suggested previously,^12,21 fits in with the first order rate dependence on the SO₃ concentration in apolar solvent. Morkovnik and Akopova²³ have proposed an alternative low-barrier mechanism with one SO₃ in apolar solvent. In their proposal, a sulfuric acid molecule acts as a catalyst by transferring the proton from benzene to SO₃via a relay-race mechanism (Scheme 2).

In this study, we aim to further elucidate the reaction mechanism and explain the experimental data. To fully account for solvation and dynamic effects, ab initio molecular dynamics (AIMD) is the preferred method of choice. The article is organized as follows. We start with the proposal of a new kinetic model, which is validated by a variety of AIMD simulations in gas phase as well as fully solvated in CCl₃F and CH₃NO₂ solvent. Metadynamics (MTD) simulations are performed to study the stability and reactivity of benzene–SO₃ and benzene–S₂O₆ σ-complexes, and to estimate the free energy surface (FES) of the sulfonation reaction. From unbiased MD simulations, we analyze intermolecular interactions between benzene, SO₃, and solvent, and we investigate the stability of benzene⋯SO₃ and benzene⋯SO₃⋯SO₃ π-complexes in each environment. Finally, restrained MD (rMD) simulations are performed to evaluate the reactivity of the benzene–S₂O₆ σ-complex, and of benzene–SO₃ + 1 catalytic H₂SO₄.

Results and discussion

Mechanism and kinetic model

To enhance readability of the study, we start with the proposal of a new sequence of reaction steps given in Scheme 3 based on our calculation results, which are presented in the following sections, and in agreement with experimental data.

Scheme 3 New proposed kinetic model for the primary sulfonation stage with SO₃.

Here, step 1 represents formation of ArH⋯SO₃, a π-complex between the arene and a first SO₃ molecule (hereafter named the primary SO₃), step 2 is the formation of the ArH⋯SO₃⋯SO₃ π-complex between ArH⋯SO₃ and a second SO₃ (the assisting SO₃), and step 3 is the actual sulfonation reaction with formation of arenepyrosulfonic acid ArS₂O₆H, the main reaction product. Step 3 requires a significant reorientation of the two SO₃ molecules toward the plane of the benzene ring to make a cyclic proton rearrangement possible. This mechanism differs from the proposed mechanism of Cerfontain and coworkers^13–16 in that it does not account for an intermediate σ-complex state. Although its presence is not ruled out, our simulations confirm previous experimental and theoretical data^12,24 that the stability of σ-complexes is generally low and may be discarded from the kinetic model. Using steady-state conditions, the reaction rate is given by the rate equation (see Appendix):

Depending on the solvent, the reaction may be first order or second order in SO₃,^13,14 which dictates limiting condition 1 (Table 1). In apolar CCl₃F solvent, the reaction is first order in SO₃, and thus k₋₁ ≪ k₂[SO₃], which implies that the ArH⋯SO₃ complex must remain stable long enough that it can associate with an assisting SO₃ before the primary SO₃ dissociates from the arene. In polar CH₃NO₂, the reaction is second order in SO₃ and k₋₁ ≫ k₂[SO₃], thus the arene exists as mostly uncomplexed with SO₃. Whereas the rate-limiting step could in principle be step 1 or 3 in CCl₃F, and step 2 or 3 in CH₃NO₂ (Table 1, condition 2), our results suggest that (i) the calculated free energy barriers of sulfonation are relatively high, and (ii) the assisting SO₃ detaches rapidly from ArH⋯SO₃⋯SO₃, which strongly suggests that k₋₂ ≫ k₃, thus step 3 is rate-limiting in both solvents.

Table 1 Experimental order in SO₃ (x) in solvents CCl₃F and CH₃NO₂, corresponding limiting conditions and rate-limiting step (RLS) derived from eqn (1)

Solvent	[SO3]^x	Condition 1	Condition 2	RLS
CCl3F	1	k −1 ≪ k2[SO3]	k −2 ≪ k3	1
			k −2 ≫ k3	3
CH3NO2	2	k −1 ≫ k2[SO3]	k −2 ≪ k3	2
			k −2 ≫ k3	3

---

In the following sections we validate the proposed mechanism and kinetic model by analyzing a carefully chosen set of advanced AIMD simulations in gas phase and fully solvated in CCl₃F and CH₃NO₂.

Reaction of benzene + 1 SO₃

We first examine if a stable σ-complex can be formed between benzene and a single SO₃ molecule, as proposed by Cerfontain and coworkers,^13–16 and if aromatic sulfonation is possible with only 1 SO₃. MTD simulations are performed with one collective variable (CV), the coordination number CN_CS that measures bond formation between a benzene carbon and the SO₃ sulfur. The resulting potential of mean force (PMF) in gas phase, in CCl₃F, and in CH₃NO₂ is shown in Fig. 1. The benzene⋯SO₃ π-complex is located at CN_CS ≅ 0.085 (r_CS = 2.9 Å). Large free energy differences are observed between the three environments, the σ-complex (CN_CS ≅ 0.5–0.6) being most stabilized in the polar CH₃NO₂. However, no intermediate σ-complex between benzene and SO₃ can be discerned from the PMF, in agreement with previous static calculations.¹²

Fig. 1 Reaction between benzene and 1 SO₃ in gas phase, CCl₃F, and CH₃NO₂. PMF of 1D MTD with error bars representing the standard deviation. Numbers above the curves indicate free energy ΔF (kJ mol⁻¹) of the σ-complex at CN_CS = 0.6 relative to the reactant state. Structures refer to the gas phase curve.

The MTD simulations are continued by further adding Gaussian hills to the reactant state until reaction. No reaction is observed until a barrier height ΔF^‡ ≅ 244 kJ mol⁻¹ is reached, which is incompatible with the fast experimental reaction kinetics.¹⁵ In conclusion, the MTD simulations suggest that sulfonation of benzene with a single SO₃ is not a viable pathway, and can be discarded in the kinetic analysis.

Stability of benzene⋯SO₃⋯SO₃ and benzene⋯SO₃ π-complexes

To assess the validity of steps 1 and 2 of our kinetic model (Scheme 3), the stability of the π-complexes is studied with a series of independent (different starting structures and velocities) and unbiased AIMD simulations starting from the benzene⋯SO₃⋯SO₃ complex. The results are listed in Table 2 and summarized in Fig. 2.

Table 2 Overview of MD and rMD simulations

Initial conditions	Environment	t (ps)	Reaction^b
a Total simulation time. b Time at which the event takes place is given between parentheses.
Unbiased MD
Benzene⋯SO3⋯SO3	Gas phase	120	Remains stable
	CCl3F	200	→Benzene⋯SO3 + SO3 (30 ps)
		200	→Benzene⋯SO3 + SO3 (140 ps)
	CH3NO2	70	→Benzene⋯SO3 + SO3 (30 ps) →Benzene + 2 SO3 (70 ps)
		180	→Benzene⋯SO3 + SO3 (60 ps) →Benzene + 2 SO3 (180 ps)
Restrained MD
Benzene–S2O6	Gas phase	40	→Benzene–S2O6H
		60	→Benzene–S2O6H
	CCl3F	5	→Benzene–S2O6H
		10	→Benzene–S2O6H
	CH3NO2	200	Remains stable
Benzene–SO3⋯H2SO4	Gas phase	6	→Benzene–SO3H⋯H2SO4
	CCl3F	170	→Benzene–SO3H⋯H2SO4
	CH3NO2	90	→Benzene–SO3H⋯H2SO4

---

Fig. 2 Stability of π-complexes in different solvation models, showing the decrease in number of SO₃ molecules complexed with benzene.

In gas phase, the benzene⋯SO₃⋯SO₃ complex remains stable throughout the simulation. A clear distinction can be made between the two interacting SO₃ molecules (Fig. 3a). The first (primary) SO₃ is tightly bound to benzene by strong π interaction between its S atom and the π-conjugated ring system located above the benzene C atoms. The second (assisting) SO₃ is loosely bound to the primary SO₃ and occasionally strays away from the complex, while remaining at the same side of the ring plane most of the time. In both CCl₃F and CH₃NO₂ solvents, however, one SO₃ dissociates from the complex and drifts into the bulk solvent. The resulting benzene⋯SO₃ complex remains stable in CCl₃F, whereas in CH₃NO₂ this complex dissociates into solvated benzene and SO₃.

Fig. 3 Intermolecular interactions observed during AIMD simulations. (a) Benzene–SO₃ interactions in gas phase. Superposition of snapshots of the benzene⋯SO₃⋯SO₃ complex in top view (left) and side view (right). Dots represent the S atom positions of the primary SO₃ (red) and assisting SO₃ (blue). (b) Specific SO₃–CH₃NO₂ interactions. Probability isosurface (red) of CH₃NO₂ O atoms around the SO₃ S atom in a benzene⋯SO₃ complex (left) and in separated SO₃ (right). (c) Specific benzene–CH₃NO₂ interactions. Probability isosurface of solvent O (red) and C (grey) atoms around benzene C atoms, in a benzene⋯SO₃ complex (left) and in separated benzene (right). Deviations from symmetry in the isosurfaces are due to finite sampling.

Overall, the unbiased MD simulations described here show that the benzene⋯SO₃⋯SO₃ π-complex is stable in gas phase but unstable in both solvents. The benzene⋯SO₃ complex appears relatively stable in CCl₃F, consistent with the limiting condition k₋₁ ≪ k₂[SO₃] (Table 1), and in agreement with the experimental first order rate in SO₃. In contrast, strong interactions with CH₃NO₂ solvent shift the equilibrium toward separate benzene and SO₃ reactants, in accord with k₋₁ ≫ k₂[SO₃] (Table 1), in correspondence with the second order rate in SO₃. In summary, our simulation results suggest that the reactant species is benzene⋯SO₃ in CCl₃F and benzene CH₃NO₂, in agreement with steps 1 and 2 of our kinetic model (Scheme 3).

Specific solvation of reactants

The low stability of the benzene⋯SO₃ π-complex in CH₃NO₂ is attributed to specific solvent interactions. Electrostatic interaction between the positively charged S atom of SO₃ and the negatively charged O atoms of CH₃NO₂ strongly increases upon dissociation of the complex. This solvent interaction is demonstrated in Fig. 3b by the red probability isosurface of CH₃NO₂ O atoms around SO₃, where larger volumes indicate a higher probability of finding an O atom inside the isosurface. Likewise, electrostatic interactions between the benzene π electron cloud and the positively charged methyl group of CH₃NO₂ (grey), and between the positively charged H atoms of benzene and the O atoms of CH₃NO₂ (red) markedly increase upon dissociation (Fig. 3c).

Reaction of benzene + 2 SO₃

The sulfonation reaction starting from the benzene⋯SO₃⋯SO₃ π-complex, corresponding to step 3 in our kinetic model (Scheme 3), is investigated with MTD and two collective variables, CN_CS and CN_CO; the latter is the coordination between a benzene C and the O atom of the assisting SO₃. In Fig. 4a–c, the resulting 2D FES is shown as the average of three MTD simulations. Again, no intermediate σ-complex is observed in the three environments. The surface is very similar in gas phase and CCl₃F, whereas in CH₃NO₂ the reactant valley is much shallower due to stabilization of TS by the electrostatic field. The position of TS, as indicated with white dots, is similar in the three environments: CV1 ≅ CV2 ≅ 0.6, which corresponds to a cyclic transition state with two SO₃ molecules at r_CS ≅ 1.84 Å and r_CO ≅ 2.7–3.2 Å. In all cases benzenepyrosulfonic acid is formed in one step. In contrast to CV1, the free energy change along the CV2 axis is small. After projecting CV2 onto CV1, a 1D PMF is obtained (Fig. 4d). Mean geometric parameters at TS are shown in Table S2 (ESI†).

Fig. 4 FES of 2D MTD in gas phase (a), CCl₃F (b), and CH₃NO₂ (c), and corresponding projected PMFs along CN_CS (d). The given numbers indicate average ΔF^‡ values (kJ mol⁻¹).

Low-temperature MD of the σ-complex in CH₃NO₂ shows that a metastable intermediate state may exist nonetheless (see ESI†). The stability of the σ-complex is however very low, which may in part be due to the BLYP density functional used in this study. Low stabilities of the intermediate state (2–4 kJ mol⁻¹) have also been calculated previously at the M06-2X/6-311+G(2d,2p) level for benzene and 1,4-dichlorobenzene in implicit CH₃NO₂.¹²

Reactivity of benzene–S₂O₆ σ-complex

Although we have now established that σ-complexes of benzene and SO₃ are fairly unstable, it is insightful to analyze the benzene–SO₃ interactions in the σ-complex state and to monitor the influence of the environment on their reactivity toward benzene–S₂O₆H formation. Additionally, the obtained reactivity information will allow us to decompose the free energy barrier into several discrete contributions.

Restrained MD (rMD) simulations are performed of the benzene–SO₃ σ-complex + an assisting SO₃, which forms a benzene–S₂O₆ σ-complex. To stabilize the σ-complex, the C₁–S₁ bond is restrained with a harmonic potential, which keeps CN_CS between 0.5 and 0.6 until sulfonation takes place. The results are summarized in Table 2. In gas phase (40–60 ps) and in CCl₃F solvent (5–10 ps), a concerted sulfonation reaction takes place by proton transfer to the assisting SO₃ with formation of benzenepyrosulfonic acid in one step (Fig. 5). In CH₃NO₂ solvent, however, no product formation is observed during 200 ps. Two factors contribute to the increased stability of the σ-complex in CH₃NO₂: (i) stabilization of the zwitterionic TS by the strong solvent electrostatic field, and (ii) competition with the solvent for hydrogen bonding with H₁. The latter effect is demonstrated in Fig. 6a, showing a stable intramolecular O₂–H₁ hydrogen bond in gas phase and CCl₃F, whereas in CH₃NO₂ this hydrogen bond is present only ∼50% of the time.

Fig. 5 TS of sulfonation of benzene with 2 SO₃. In CCl₃F, the reaction takes place spontaneously during rMD. In CH₃NO₂ an additional barrier needs be crossed with rMTD. Atom colors are: C (grey), Cl (green), F (cyan), H (white), N (blue), O (red), S (yellow).

Fig. 6 Geometric probability densities (P) in reactant state (benzene⋯SO₃ π-complex, full lines) and close to TS (benzene–S₂O₆ σ-complex, dashed lines), in gas phase (blue), CCl₃F (green), and CH₃NO₂ (red). (a) minimal distance r(H₁–O₂) and minimal distance r(H₁–O) with CH₃NO₂ (filled grey), (b) angle α(C₄–C₁–H₁), (c) HOMA aromaticity index, (d) distance r(C₁–H₁).

Solvation effects at the transition state

The segments of rMD trajectories before reaction are analyzed to determine how the transition state geometry is affected by solvation effects. The effects of the solvent on benzene geometry near TS are striking (Fig. 6b–d). The angle α(C₄–C₁–H₁) is calculated as a measure for the degree of sp³ hybridization at C₁ and thus the stability of the σ-complex. In the reactant state, the mean angle 〈α〉 = 174° in all three environments, indicating that H₁ is located mostly in the plane of the benzene ring and hence C₁ is almost purely sp² hybridized, whereas near TS C₁ displays strong sp³ character. The smallest 〈α〉 is measured in CH₃NO₂ (131°), followed by CCl₃F (138°) and gas phase (147°). The strong stabilizing effect of CCl₃F is remarkable given its low polarity. The harmonic oscillator model of aromaticity index (HOMA)²⁵ fluctuates around a mean value 〈HOMA〉 = 0.59–0.63 in the reactant state, but degrades near TS, indicating nearly complete loss of aromaticity. This effect is especially strong in CH₃NO₂ (〈HOMA〉 = −0.19), in comparison with CCl₃F (0.02) and gas phase (0.13). The C₁–H₁ bond is stretched from 〈r_CH〉 = 1.102–1.104 Å to 1.124, 1.121, and 1.113 Å in CH₃NO₂, CCl₃F, and gas phase, respectively. In CCl₃F, a bimodal distribution is observed due to frequent near-proton transfer events taking place.

An alternative pathway in CH₃NO₂

An estimate of the additional free energy barrier required for sulfonation of the benzene-S₂O₆ σ-complex in CH₃NO₂ is obtained with restrained MTD (rMTD) (for computational details and PMFs, see ESI†), yielding an average additional ΔF = 12 kJ mol⁻¹. In one of the rMTD runs, an alternative pathway is observed (Fig. 7). First, the C–H bond is weakened by hydrogen bonding interaction with the assisting SO₃, which strongly lowers the activation free energy of proton transfer. Next, the proton is transferred to a nitromethane molecule (224 fs, TS1). Shortly thereafter, the proton is transferred from nitromethane to the primary SO₃ to form benzenesulfonic acid (404 fs, TS2). Thus, although in this mechanism the assisting SO₃ is not the proton acceptor, it still plays an essential role in the reaction by activating the C₁–H₁ bond. The calculated additional free energy barrier for this pathway ΔF^‡ = 14 kJ mol⁻¹, suggesting that this new mechanism is feasible and may compete with the cyclic proton transfer pathway.

Fig. 7 Snapshots of the new CH₃NO₂-mediated sulfonation pathway. (0 fs) H₁–O₂ hydrogen bond at maximal strength. (224 fs) TS1: proton transfer from benzene to CH₃NO₂. (404 fs) TS2: proton transfer from CH₃NO₂ to the primary SO₃. (2404 fs) stable benzenesulfonic acid product. Blue arrows indicate the direction of motion of the proton.

Free energy decomposition

Combining the results of 1D and 2D MTD and the rMD simulations, the total ΔF^‡ can be decomposed into three contributions: (1) benzene⋯SO₃ π-complex to benzene–SO₃ σ-complex, (2) benzene–SO₃ σ-complex to benzene–S₂O₆ σ-complex, and (3) benzene–S₂O₆ to TS. In gas phase, the presence of an assisting SO₃ strongly reduces the free energy required to form a σ-complex from the π-complex (177 to 129 kJ mol⁻¹). Once the benzene–S₂O₆ σ-complex is formed, sulfonation proceeds spontaneously. In CCl₃F, the assisting SO₃ has a slight stabilizing effect on the σ-complex (129 to 125 kJ mol⁻¹). In CH₃NO₂, the solvent already has a strong stabilizing effect, which causes the additional stabilizing effect of the assisting SO₃ to be rather small (85 to 76 kJ mol⁻¹). The combined stabilizing effects lead to an additional free energy barrier required for sulfonation (76 + 12 = 88 kJ mol⁻¹).

H₂SO₄ as a catalyst

As an alternative explanation for the first-order reaction in SO₃ in noncomplexing media, Morkovnik and Akopova suggested a relay-race mechanism involving a brønsted acid as a catalyst (Scheme 2).²³ The plausibility of this mechanism is verified with rMD simulations of the benzene–SO₃ σ-complex, restrained at the C₁–S₁ bond, in the presence of a H₂SO₄ molecule. Spontaneous sulfonation takes place in all three environments within 170 ps (Table 2), although the mechanism differs from the mechanism proposed previously²³ and depends on the polarity of the environment (Fig. 8). In gas phase and CCl₃F, H₁ is first transferred from benzene to H₂SO₄ at TS1, stabilized by one (gas phase) or two (CCl₃F) H₂SO₄–SO₃ hydrogen bonds. Proton transfer from H₂SO₄ to SO₃ follows quickly thereafter (TS2, 16 fs). In CH₃NO₂, first proton transfer occurs from H₂SO₄ to SO₃ (TS1), followed by benzene-to-H₂SO₄ proton transfer 106 fs later (TS2). We conclude that catalytic amounts of H₂SO₄ in CCl₃F may provide an alternative pathway consistent with first-order kinetics in SO₃, however the participation of H₂SO₄ is not required in our kinetic model. In CH₃NO₂, the stability of the benzene⋯SO₃⋯H₂SO₄ π-complex is probably too low to meaningfully contribute to the overall reaction.

Fig. 8 TS1 (top) and TS2 (bottom) in the presence of catalytic H₂SO₄, corresponding to the first and second proton transfer in CCl₃F (left) and CH₃NO₂ (right).

Conclusions

The mechanism and kinetics of electrophilic aromatic sulfonation have been investigated in great detail with various DFT-based first-principles MD and MTD simulations. Three different environments were compared: gas phase, and fully solvated in explicit apolar (CCl₃F) and polar (CH₃NO₂) solvents. Several alternative reaction mechanisms were evaluated, including with 1 or 2 SO₃ molecules, with a catalytic H₂SO₄ molecule, and with a participating solvent molecule. The kinetic model proposed by Cerfontain and coworkers, involving stable ArH–SO₃ and ArH–S₂O₆ σ-complexes in step 1 and 2, was examined. Our results suggest that both benzene–SO₃ and benzene–S₂O₆ σ-complexes are unstable at room temperature, thus ruling out the Cerfontain model. Our simulation data confirm the static gas phase and implicit solvent calculations performed by Schleyer and coworkers, which suggest that the free energy barrier for sulfonation with a single SO₃ is too high to be feasible. In all three environments, a low-energy concerted pathway starting from the benzene⋯SO₃⋯SO₃ π-complex was found, involving a cyclic transition state with proton transfer from benzene to the assisting SO₃, and with formation of benzenepyrosulfonic acid. This mechanism is in agreement with the calculations of Schleyer and coworkers, although in CH₃NO₂ an intermediate cyclic σ-complex state with low stability was found by them.

In order to bring the concerted mechanism into agreement with experimental kinetic data, a new kinetic model was proposed. Steps 1 and 2 represent the formation of ArH⋯SO₃ and ArH⋯SO₃⋯SO₃ π-complexes, respectively, followed by cyclic proton transfer with formation of benzenepyrosulfonic acid in step 3. A rate equation was calculated using the steady-state approximation, and the limiting conditions were determined on the basis of the experimental rate order in SO₃. In CCl₃F we found that k₋₁ ≪ k₂[SO₃], whereas in CH₃NO₂ the opposite condition k₋₁ ≪ k₂[SO₃] applies. The validity of steps 1 and 2 was confirmed with long timescale MD simulations, suggesting that the benzene⋯SO₃ π-complex is relatively stable in CCl₃F but quickly dissociates in CH₃NO₂. The mechanism suggested by Morkovnik and Akopova involving one SO₃ and a catalytic H₂SO₄ molecule was shown to provide a viable alternative route to sulfonation in CCl₃F. Furthermore, a new and unanticipated mechanism was discovered in CH₃NO₂, in which two SO₃ and one CH₃NO₂ cooperate in a stepwise proton transfer from benzene to CH₃NO₂, followed by proton transfer to the primary SO₃, leading directly to benzenesulfonic acid. Future research will focus on sulfonation of substituted benzenes and benzene derivatives with SO₃ to investigate whether our kinetic model is applicable to a wide range of arenes.

Computational methods

Molecular dynamics

Born–Oppenheimer molecular dynamics simulations are performed at the DFT level with the gradient-corrected BLYP functional,^26,27 the DZVP-GTH basis set,²⁸ and Grimme D3 dispersion corrections.²⁹ The BLYP functional has been shown to produce barriers for SN₂ and E₂ reactions in good agreement with high-level (MP2) calculations.³⁰ The integration time step is set at 1 fs, with snapshots taken every 2 fs. Simulations in solvent are performed in the NVT ensemble, using the canonical sampling through velocity rescaling thermostat³¹ with a time constant of 50 fs. The reactants are placed inside a periodic cubic box filled with 25 CCl₃F or 40 CH₃NO₂ molecules. The densities are set to the experimental densities of CCl₃F (1.48 g ml⁻¹) and CH₃NO₂ (1.13 g ml⁻¹) at 298 K, with corresponding box lengths 15.847/15.911 Å (CCl₃F) and 15.482/15.549 Å (CH₃NO₂), depending on the number of SO₃ molecules present. All AIMD simulations are performed with the CP2K simulation package (version 2.6).³²

Metadynamics

Most chemical reactions are not accessible in the timescale (typically < 1 ns) that can be reached with AIMD. In order to sample the relevant TS regions, enhanced sampling techniques need to be used. Metadynamics is a nonequilibrium MD method introduced by Laio and Parrinello.^33,34 In recent studies, we have successfully used the MTD technique to estimate free energy surfaces in a variety of systems, including heterogeneous catalysis and reactions in solution.^35–37 Sampling is advanced by adding Gaussian potential hills along a limited number of carefully chosen collected variables during the simulation, effectively flattening the FES by filling low-energy regions. The FES is then calculated as the opposite of the summation of the Gaussian hills.

The hills are added every 25 steps along one or two collective variables (CVs), defined by coordination numbers CN:

where the sum runs over two nonoverlapping sets of atoms i and j, r_ij is the distance between atoms i and j, and r₀ is a reference distance. CV1 is described by CN_CS between a benzene C₁ and S₁ (see Scheme 2 for atom numbering), with r₀ = 1.964 Å. In the 2D MTD simulations, CV2 corresponds to CN_CO, the sum of coordination numbers between C₁ and the three O₂ atoms, with r₀ = 2.719 Å. The width of the hills is set to 0.02. In the 1D MTD simulations, the hill height is initially set to H = 1 kJ mol⁻¹, and reduced to 0.5 kJ mol⁻¹ after 40 ps. In the 2D MTD, H = 2 kJ mol⁻¹. To limit sampling to a region close to the bound state, half harmonic bias potentials are added to CV1 and CV2 at position 0.03 with a force constant K_f = 100 a.u. In the 2D MTD, an additional half harmonic potential is added to the S₁–S₂ distance at 4.5 Å with K_f = 0.19 a.u. The MTD simulations are ended once the product state is reached. To dampen excessive proton fluctuations, the mass of H₁ is increased to that of tritium.

In the 2D MTD, Gaussian hills are placed along two CVs simultaneously: CN_CS and CN_CO, which allows for a sulfonation with either one or two SO₃ molecules. Due to the involvement of proton transfer, however, it is impossible to simulate the reaction in the conventional way, i.e. sampling many forward and reverse reactions by filling both the reactant and product wells with Gaussian hills. Instead, three independent MTD simulations are initiated from different starting geometries and terminated as soon as the proton transfer has taken place. The TS is then taken as the last stationary point along the r_CH trajectory before proton transfer takes place.

Analysis

Probability isosurfaces are calculated with the volmap tool of VMD.³⁸ The probability isosurface of CH₃NO₂ O atoms within 4 Å of the SO₃ S atom is calculated with isovalue = 65%. Probability isosurfaces of CH₃NO₂ O and C atoms within 4 Å of any of the six benzene C atoms are calculated with isovalue = 30%. Conversion of the 2D FES into a 1D PMF is achieved by projecting CV2 onto CV1, and free energy barriers are calculated as the free energy difference between the TS and the reactant state.³⁶

Appendix

Here we derive the rate equation (eqn (1)). Let A = ArH, B = SO₃, C = ArH⋯SO₃, D = ArH⋯SO₃⋯SO₃, E = ArS₂O₆H. Using the steady-state approximation, we have

From eqn (3)–(4):

The form of this equation differs from the rate equation presented by Lammertsma and Cerfontain,¹⁷ who in the denominator neglect k₋₂ in comparison to k₃ in the term (k₋₂ + k₃)k₂[B], but not in (k₋₂ + k₃)k₋₁. No rationale was given for the partial neglect of k₋₂, and it does not seem to simplify the kinetic analysis.

Acknowledgements

The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Hercules Foundation and the Flemish Government – department EWI. P. G. and F. D. P. wish to thank the Fund for Scientific Research – Flanders (FWO) and the Vrije Universiteit Brussel (VUB) for their continuous support. They also specifically want to mention the Strategic Research Program awarded to the ALGC group by the VUB which started on January 1, 2013.

References

1. E. A. Knaggs and M. J. Nepras, in Kirk-Othmer Encyclopedia of Chemical Technology, John Wiley & Sons, Inc., 2000,  DOI:10.1002/0471238961.1921120611140107.a01.
2. M. B. Smith and J. March, in March's Advanced Organic Chemistry, John Wiley & Sons, Inc., 2006,  DOI:10.1002/9780470084960.ch11, pp. 657–751.
3. J. S. Carey, D. Laffan, C. Thomson and M. T. Williams, Org. Biomol. Chem., 2006, 4, 2337–2347.
4. G. P. Dado, E. A. Knaggs and M. J. Nepras, in Kirk-Othmer Encyclopedia of Chemical Technology, John Wiley & Sons, Inc., 2000,  DOI:10.1002/0471238961.1921120611140107.a01.pub2.
5. D. W. Roberts, Org. Process Res. Dev., 2003, 7, 172–184.
6. Y. Chen, Y. Su, F. Jiao and G. Chen, RSC Adv., 2012, 2, 5637–5644.
7. C. M. Suter and A. W. Weston, in Organic Reactions, John Wiley & Sons, Inc., 2004, vol. 3, pp. 141–197.
8. B. Galabov, D. Nalbantova, P. v. R. Schleyer and H. F. Schaefer, Acc. Chem. Res., 2016, 49, 1191–1199.
9. J. Kong, B. Galabov, G. Koleva, J.-J. Zou, H. F. Schaefer and P. v. R. Schleyer, Angew. Chem., Int. Ed., 2011, 50, 6809–6813.
10. B. Galabov, G. Koleva, S. Simova, B. Hadjieva, H. F. Schaefer and P. v. R. Schleyer, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 10067–10072.
11. G. Koleva, B. Galabov, B. Hadjieva, H. F. Schaefer and P. v. R. Schleyer, Angew. Chem., Int. Ed., 2015, 54, 14123–14127.
12. G. Koleva, B. Galabov, J. Kong, H. F. Schaefer and P. v. R. Schleyer, J. Am. Chem. Soc., 2011, 133, 19094–19101.
13. J. K. Bosscher and H. Cerfontain, Recl. Trav. Chim. Pays-Bas, 2010, 87, 873–887.
14. J. K. Bosscher and H. Cerfontain, Tetrahedron, 1968, 24, 6543–6555.
15. J. K. Bosscher and H. Cerfontain, J. Chem. Soc. B, 1968, 1524–1526,  10.1039/j29680001524.
16. H. Cerfontain, Recl. Trav. Chim. Pays-Bas, 2010, 104, 153–165.
17. K. Lammertsma and H. Cerfontain, J. Chem. Soc., Perkin Trans. 2, 1980, 28–32,  10.1039/p29800000028.
18. H. Cerfontain and A. Koeberg-Telder, Recl. Trav. Chim. Pays-Bas, 2010, 89, 569–574.
19. A. Koeberg-Telder and H. Cerfontain, Recl. Trav. Chim. Pays-Bas, 2010, 90, 193–206.
20. J. O. Morley and D. W. Roberts, J. Org. Chem., 1997, 62, 7358–7363.
21. J. O. Morley, D. W. Roberts and S. P. Watson, J. Chem. Soc., Perkin Trans. 2, 2002, 538–544,  10.1039/b109338j.
22. F. H. Westheimer, Chem. Rev., 1961, 61, 265–273.
23. A. S. Morkovnik and A. R. Akopova, Dokl. Chem., 2013, 450, 122–126.
24. F. Kučera and J. Jančář, Polym. Eng. Sci., 1998, 38, 783–792.
25. T. M. Krygowski and M. K. Cyranski, Phys. Chem. Chem. Phys., 2004, 6, 249–255.
26. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
27. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
28. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter, 1996, 54, 1703–1710.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
30. B. Ensing and M. L. Klein, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6755–6759.
31. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
32. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
33. A. Laio and F. L. Gervasio, Rep. Prog. Phys., 2008, 71, 126601.
34. A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
35. S. L. C. Moors, B. Brigou, D. Hertsen, B. Pinter, P. Geerlings, V. Van Speybroeck, S. Catak and F. De Proft, J. Org. Chem., 2016, 81, 1635–1644.
36. S. L. C. Moors, K. De Wispelaere, J. Van der Mynsbrugge, M. Waroquier and V. Van Speybroeck, ACS Catal., 2013, 3, 2556–2567.
37. J. Van der Mynsbrugge, S. L. C. Moors, K. De Wispelaere and V. Van Speybroeck, ChemCatChem, 2014, 6, 1906–1918.
38. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38 , 27-38.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6sc03500k

---