Mechanism, kinetics and selectivity of a Williamson ether synthesis: elucidation under different reaction conditions

The best route to uncover the mechanism of chemical reactions remains a topic of intense debate. In this work, we deploy a three-faceted approach that combines experimental probing, detailed kinetic modelling and quantum-mechanical calculations for the study of the mechanism and regioselectivity of a Williamson ether synthesis, which is of interest because of its simplicity and its broad scope in laboratory and industrial synthesis. The choice of solvent is found to have a large impact on the experimental regioselectivity, with ratios of O -alkylated to C -alkylated product at 298 K of 97:3 in acetonitrile and of 72:28 in methanol. Through experiments and kinetic modelling, we identify reaction networks that differ significantly from solvent to solvent, providing insights into the factors (proton-exchange, solvolysis and product degradation) that impact on regioselectivity and the relative rates of O -alkylation and C -alkylation. The kinetic models yield detailed information on reaction rates and energy barriers and on the existence of an additional double alkylation pathway. We carry out quantum mechanical calculations and elucidate the transition states for the two main alkylation pathways. The quantum-mechanical calculations highlight structural differences between the transition states found for the two alkylation pathways and provide information on the effect of the solvent on the stabilisation/destabilisation of various structures and hence on reaction selectivity. The three-faceted approach provides complementary information into the elementary steps of the reaction mechanism.


Introduction
Understanding the mechanism of a chemical reaction and quantifying reaction progress as a function of conditions are crucial steps in numerous practical applications in organic chemistry, 1 for example when optimizing reaction conditions with respect to yield and selectivity. Postulating a mechanism requires accounting for the elementary steps of the reaction, the nature of intermediates and transition states, any side reactions, the energetics of transformations and the rate laws describing each reaction step (i.e., reaction kinetics). External factors, such as temperature, pressure, solvent media and catalysts, as well as internal factors, such as reactant concentrations, physical states and surface areas, can affect the reaction mechanism, making it difficult to achieve a quantitative description of the process, especially when multiple steps and alternative pathways are present. 2 In view of the complexity of reaction mechanisms, efforts have been made to combine experimental and computational tools to achieve a detailed understanding of reaction mechanisms and derive predictive results, but the role of different tools remains a subject of debate for the scientific community. Experimental techniques are most commonly employed in the study of reaction mechanisms and quantum-mechanical (QM) methods for electronic-structure calculations are often seen to be valuable in understanding of reaction mechanisms. 1 Several authors have questioned the value of using QM methods in isolation to investigate reaction mechanisms. [2][3][4] Instead, the synergistic combination of experiments and QM calculations has been suggested as an effective approach to support or disprove a proposed mechanism. 1,5,6 The benefits of combining experimental and QM computational approaches for the elucidation of reaction mechanisms have been illustrated in recent reviews. 1,[6][7][8][9][10] In particular, Plata and Singleton 2 have shown that in the case of an alcoholmediated Morita Baylis-Hillman reaction, QM computational studies carried out in isolation and focused on the simplest mechanism fail to capture the actual reaction mechanism. On the other hand, calculations carried out in concert with an extensive experimental study provide insights into steps that cannot be probed experimentally, e.g., the aldol step. Liu et al. 11 have shown that in the case of transition-metal catalyzed Negishi coupling reaction, in situ IR experiments reveal that the reaction is promoted via a transmetalation path that occurs competitively to the reductive elimination path, which is traditionally considered as the key path leading to the cross-coupling products in the Negishi coupling reaction. Computations of energy barriers and substitution effects for the two competing paths have provided corroborating evidence in favour of the transmetalation pathway, resulting in the design of an experimental strategy to avoid the homocoupling and dehalogenation side products that hinder the scale-up of this reaction.
Singleton and Wang 12 have shown that the identification of the best operative transition state (TS) for the enantioselective epoxidation of trans-β-methylstyrene is largely facilitated by the use of experimental values for the kinetic isotope effect (KIE). The initial conformational space of 66 calculated TSs was reduced to 3 TS structures through comparison of experimental and predicted KIE values. For the TS with the lowest calculated energy, the computed KIE values are in excellent agreement with the experimental values. These examples illustrate the wealth of mechanistic information that is uncovered when experiments and QM computations are employed synergistically. The combination of experimental and QM computational approaches can be particularly useful to enhance the understanding of fundamental chemical problems, 1 such as the effects of additives, bases and acids on reactivity, 13 the effects of solvents on reactivity and selectivity 14 and the effects of ligands on bond and coordination energy. [15][16][17] To gain a more complete understanding of a reaction, experimental and QM approaches can be complemented by a third means of analysis, namely kinetic modelling. Timeresolved experimental techniques enable the development of kinetic models that describe the observed reaction kinetics mathematically. Kinetic modelling adds a new dimension as it can be used to predict the evolution of species concentrations under conditions not assessed experimentally. Moreover, such models help to assess the validity of the mechanism postulated for a given reaction and provide a vital, quantitative link between experimental data and QM models. The role of kinetic modelling in reaction mechanism investigation has for instance been demonstrated through the use of reaction progress kinetic analysis (RPKA). 18 In the study of the palladium-catalyzed Heck coupling reaction of aryl halides, 18 where RPKA was applied, the combination of just three time-series experiments with different initial substrate concentrations and two graphical reaction-rate equations revealed that neither catalyst deactivation nor product inhibition influence the reaction. Such an effective analysis of the reaction mechanism was possible thanks to the abundance of experimentally measured data available to construct the graphical reaction-rate equations. Despite the usefulness of kinetic models in understanding complex reactions and in supporting subsequent process development, the derivation of such models is too often overlooked in mechanistic studies.
In this work, we explore the value of following a threefaceted approach that combines experiments, kinetic modelling and QM calculations by analysing the impact of reaction conditions on the mechanism of an ether synthesis reaction. We investigate how each method of analysis can inform and support the findings derived from the other two by focusing on a relatively simple and well-established ether synthesis, the Williamson reaction between sodium β-naphthoxide and benzyl bromide, shown in Scheme 1. We choose to focus on a Williamson ether synthesis due to the industrial and scientific relevance of this class of reactions. C-O bond formation is necessary for the manufacturing of important intermediates and products 19 and Williamson reactions are used extensively industrially. 20,21 Williamson reactions are also used routinely in pharmaceutical studies, for example in the synthesis of anti-influenza agents. 22 As they involve a bimolecular nucleophilic substitution at a saturated carbon, they are representative of one of the most widely used reaction classes in synthetic organic chemistry. The particular example under study addresses the question of regioselectivity when an ambident nucleophile is used; this is a common scenario in synthetic chemistry where the effect of solvent on product ratio can be difficult to predict using standard heuristics.
The solvent medium is known 23 to impact the kinetics and the selectivity of the reaction in forming either benzyl β-naphthyl ether (O-alkylated product) or 1-benzyl-2-naphthol (C-alkylated product). The final C-alkylated product is initially formed as the keto isomer of the product which quickly tautomerizes to the more stable phenol form (Scheme 2), subsequently converting to the final C-alkylated product through a proton exchange equilibrium (Scheme 3). Product ratios in different solvents have been reported in several experimental studies. [23][24][25][26][27][28] Moreover, under some circumstances, 23 the formation of an additional side product, 1,1-dibenzyl-2-(1H)-naphthalenone, is observed; this is a double C-alkylated species resulting from further alkylation of the C-alkylated product with benzyl bromide, shown in Scheme 4. To the best of our knowledge, no previous computational studies on this specific Williamson reaction have been carried out and the only mechanistic postulate available is based on a pictorial description of the transition states. 23 We perform a detailed investigation of the mechanism of this reaction, presenting new in situ NMR data acquired over a range of temperatures and concentrations in a polar aprotic solvent (acetonitrile) and a protic solvent (methanol), as well as a mechanistic and thermodynamic analysis based on these data. We first summarize the methods used in the investigation. We then present the derivation of kinetic models for the reaction in acetonitrile and methanol, based on extensive experimental data, highlighting differences in mechanism and energetics. The observed differences are elucidated via QM calculations of the key TS structures. The suitability of QM calculations to quantify energetic differences is also assessed.

Experimental methods
2.1.1 Spectroscopic methods. One-dimensional (1D) 13 C NMR spectroscopy is used to elucidate chemical structure differences among reactants and products and in situ 1D 1 H NMR is used to monitor the kinetics of Scheme 1. Twodimensional (2D) NMR spectroscopy, specifically 1 H-13 C heteronuclear single quantum correlation ( 1 H-13 C HSQC) NMR, is employed for structural assignment through a map of correlations of directly-bonded proton and carbon atoms, and 1 H-13 C heteronuclear multiple bond correlation ( 1 H-13 C HMBC) NMR is used to capture correlations between proton and carbon nuclei over longer distances while suppressing the direct-bond correlations. Control experiments to confirm the reversibility of reactions in Scheme 1 were performed in previous work, 28 mixing product species 3 and 4 together and showing that the two products are stable and no reaction occurred. A 500MHz Bruker CryoProbe Prodigy 29 spectrometer with a broadband cryogenic cooler probe, with an open-cycle liquid nitrogen cooling system for enhanced sensitivity, is used. In situ 1 H NMR spectroscopy is suitable for monitoring the kinetics of Scheme 1 as the two main products have distinguishable aliphatic CH 2 peaks in the 1 H NMR spectrum. We use equimolar amounts of the reactants, with an initial concentration of 0.09 M each. 1,3,5-Trimethoxybenzene is chosen as internal standard, with characteristic peaks at 3.74 ppm and 6.10 ppm that do not overlap in the 1 H NMR spectrum with any other peak of the species in Scheme 1. Prior to each experiment we weigh the amounts of the two reactants using an electronic scale. We then dissolve the sodium β-naphthoxide into the solution of internal standard and the solvent in a vial and transfer the content of the vial to an NMR tube with a syringe. We add benzyl bromide into the NMR tube to initiate the reaction and we transfer the tube into the NMR instrument, which is already tuned to the selected temperature. We monitor the reaction kinetics in deuterated acetonitrile (acetonitrile-d 3 ) at temperatures 298, 313 and 323 K, and in deuterated methanol (methanol-d 4 ) at temperatures 298, 313 and 328 K. In the experiments in acetonitrile, due to the low signal-to-noise ratio, no concentration data were collected prior to 2500 s for the C-alkylated species. Further experimental details (e.g., reagent source and purity) can be found in the ESI. † 2.1.2 Proton-exchange measurements. At 298 K the calculated aqueous pK a values of β-naphthol 30 and the C-alkylated product, 31 9.57 ± 0.10 and 9.60 ± 0.50, respectively, suggest that β-naphthol is not significantly more acidic than the C-alkylated product, so that proton exchange may occur between sodium β-naphthoxide and the C-alkylated product. Additionally, with the naphthoxide involved as the predominant base species in the tautomerisation of the keto isomer of the C-alkylated product, which results in naphthol and the deprotonated form of the C-alkylated product, the two species are also prone to proton exchange. We test this hypothesis by Scheme 2 The tautomerisation reaction of the initially formed keto tautomer of the C-alkylated product (4*) towards the more stable phenol form (5), promoted by the predominant base species present, naphthoxide.

Scheme 4
The reaction of benzyl bromide with the C-alkylated product (deprotonated) resulting in the formation of 1,1-dibenzyl-2-(1H)-naphthalenone (double C-alkylated product). using 13 C NMR spectroscopy to follow the aromatic carbon (C 2 ) in sodium β-naphthoxide, in β-naphthol and in the C-alkylated product, as shown in Fig. 1. Pure samples of each of the three components and a sample of a mixture of sodium β-naphthoxide and the C-alkylated product in a 1 : 1 molar ratio are prepared in DMSO. We anticipate different chemical shifts for the aromatic carbon in each one of the samples tested, depending on the protonation state of the atoms attached to it. Detection of this proton-exchange phenomenon using 1 H NMR would be very difficult, as the exchange happens very fast and the resulting peaks for the -OH groups in β-naphthol and the C-alkylated product are indistinguishable.

Estimation of kinetic parameters.
To investigate different reaction schemes, we formulate concentration-based kinetic models as a mixed set of ordinary differential equations and algebraic equations. Values of the kinetic parameters (reaction rate and Arrhenius constants, activation energies) are obtained by regression to the in situ experimental concentration data acquired via 1 H NMR spectrometry, using the parameter estimation facility in the gPROMS 32 software. The maximum likelihood formulation 33 is used with a constant variance, ranging from 0.0001 to 0.01 mol dm −3 , which are representative error values for the various concentration measurements. For the tautomerisation reaction, we apply a pseudo-steady state hypothesis based on Scheme 2 being fast so that the concentration of species 4* is kept constant. The obtained values for the corresponding estimated rate constant in the various models developed are of the order of 1000 dm 3 mol −1 s −1 and are reported in the ESI. † To assess the quality of the models, the mean average percentage error (MAPE) in concentration is reported. It is given by: where N is the number of species for which concentration data are available, N j is the number of concentration points available for species j, C i,j is the calculated concentration of species j at time point i and C exp i,j is the corresponding measured concentration. Also, the root-mean-square error (RMSE) is reported, given by: Additionally, we compare the calculated and experimental selectivities based on the percentage rate-constant ratio (PRCR). In the case where the only two products are the O-alkylated and C-alkylated products, PRCR is calculated as: where k O is the reaction rate constant of the O-alkylation reaction in Scheme 1, k C is the reaction rate constant of the C-alkylation reaction in Scheme 1 and "percentage" refers to the use of a basis of 100 for the sum of the two elements of the ratio, enabling easier comparison.
When the double C-alkylated product is present, PRCR is given by: where k DC is the reaction rate constant for the formation of double C-alkylated product. Furthermore, the computed PRCR is compared to an experimental PRCR obtained by using the measured (final) concentrations of the corresponding species, given by: where [O] and [C] are the concentration values of the O-alkylated and the C-alkylated products. When the double C-alkylated product is present, it is given by: where [DC] is the concentration value of the double C-alkylated product. 2.2.2 Density functional theory calculations. All electronicstructure calculations are performed using the software Gaussian 09 (ref. 34) (release C.01). For the reactions in Scheme 1, geometry optimizations of reactants and products and transition state searches are carried out using four density functional theory (DFT) methods (B3LYP, M05-2X, M06-2X, wB97xD) and eight basis sets (6-31G(d), 6-31+G(d), 6-31+G(d,p), 6-311+G(d,p), 6-311++G(d,p), 6-311G(2d,d,p), cc-pVDZ, cc-pVTZ) in various combinations. Frequency calculations are performed for all reaction species at the same level of theory at which they have been optimized. Vibrational frequency analyses are carried out to confirm that all reactants and products have zero imaginary frequencies and all transition states are first-order saddle points with one imaginary frequency. Gas-phase partition functions are used in the reaction rate constant calculations. The keywords ultrafine (for the integration grid), vtight (for convergence View Article Online criteria) and calcall (for analytical second derivatives at every optimization step) are used. No scaling factors are considered when calculating the partition functions and zero-point energies, except when using the B3LYP/6-31+G(d) level of theory, where a scale factor of 0.9614 is applied. 35 For the O-alkylation pathway, among the two transition-state structures identified, TS O1 and TS O2 , we select TS O2 as the one with the lowest electronic energy, assuming that the reaction proceeds only via the most favorable conformation.
To account for the presence of solvent, the continuum solvation SMD model 36 is employed, in which the solute geometry is optimized in the field generated by the dielectric constant of the solvent. The default property values of the SMD model are used for calculations in acetonitrile and methanol. The standard-state molar free energy of solvation ΔG _°, solv is calculated as: where calculations in the gas phase (denoted by IG) and liquid phase (denoted by L) are performed using optimized geometries in their corresponding phases. The term ΔE _ el,L is an electrostatic term that accounts for the change in the electronic energy when transferring a solute from the gas to the liquid phase, and the term G _ CDS,L is a non-electrostatic term which accounts for free-energy changes due to shortrange electrostatic effects and non-electrostatic effects, such as cavitation (C), dispersion (D) and solvent-related (S) structural changes. The SMD model is known to predict energies of solvation well, although larger deviations can occur when the level of theory or the solvent employed for the calculation were not included in the parametrization. 36

Reaction rate-constant calculations. For a bimolecular reaction, given as
we calculate the reaction rate constants in the liquid phase based on the transition-state theory [37][38][39] and the SMD model, 36 as where κ is the transmission coefficient based on the Wigner scheme 40 to correct for tunnelling effects, k B is the Boltzmann constant, T is the absolute temperature, h is the Planck constant, R is the ideal gas constant, c°; L i is the standard-state molar concentration for species i = A, B (reactants) and AB ‡ (transition state) for eqn (8), q ′ ;IG i T ð Þ is the ideal-gas molecular partition function of species i at temperature T (excluding the electronic energy), v i is the stoichiometric coefficient of species i with a value of −1 for each of the reactants and +1 for the transition state, ΔG _ p→c,IG i accounts for the conversion from the gas-phase standard state of T = 298.15 K and p°= 1 atm to T = 298.15 K and c°; IG i = 1 mol dm −3 as recommended by Ho et al., 41 E _ el,IG is the ideal-gas electronic energy, and Δ ‡ stands for the difference between the transition state and the reactants, for example For competing reactions with two possible pathways, the selectivity ratio of one pathway (I) versus the other (II) may be determined by the ratio of the corresponding reaction rate constants, as Assuming that the two different transition states stem from the same reactant configuration, eqn (11) may be further simplified, resulting in where ΔΔ ‡ G̲°; L II−I corresponds to the standard-state Gibbs freeenergy difference between the transition-state structures of the two alternative pathways, given as 3 Results and discussion

Reaction networks considered
We investigate three reaction networks constructed by considering different side reactions. All the reaction networks include the two reactions between sodium β-naphthoxide and benzyl bromide towards the formation of the O-and C-alkylated products (Scheme 1), the tautomerisation reaction towards the more stable form of the C-alkylated product (Scheme 2), and the equilibrium between protonated and deprotonated forms of sodium β-naphthoxide and the C-alkylated product (Scheme 3). In addition, some networks explore the importance of two further reactions: a reaction between benzyl bromide and the C-alkylated product (Scheme 4) and a reaction between benzyl bromide and methanol-d 4 (Scheme 5). The specific combinations of reactions and equilibrium we consider are presented in Table 1, with each combination corresponding to a different model. The differential equations that constitute each model are given in the ESI. † Additionally to the models presented in the main body of this study, a so-called minimal model is also developed based on the simplified Scheme shown in the ESI, † assuming that the two main O-and C-alkylated products are derived directly from the starting components while the keto isomer of the C-alkylated product is omitted. The minimal model is mentioned in the rest of the study mainly for comparison purposes.

Occurrence of proton exchange
Our study commences with an experimental investigation of the occurrence of proton exchange between sodium β-naphthoxide and the C-alkylated product, a phenomenon which has not been explicitly considered in previous studies, [23][24][25][26][27] and which we postulate on the basis of the pK a values 30 of the reaction species (i.e., β-naphthol and the C-alkylated product), as well as the resulting species of Scheme 2. The chemical shifts obtained from 13 C NMR are presented in Table 2. The displacement of the chemical shift of C 2 ( Fig. 1) in the mixture sample indicates that the magnetic shielding of the two carbon nuclei varies due to changes in the environment of the neighbouring oxygen atom in sodium β-naphthoxide and the C-alkylated product. The resulting chemical shift value for C 2 in sodium β-naphthoxide in the mixture (163.6 ppm) is approximately the average of chemical shift values for C 2 in the pure samples of naphthol and sodium β-naphthoxide. These findings validate the hypothesis that protonated and deprotonated forms of species 1 and 4 coexist in equilibrium, as per Scheme 3. We will evaluate the impact of this new finding on the time-evolution of species when performing kinetic analysis in different solvents.

Kinetics and reaction network in acetonitrile-d 3
We monitor the kinetics of the reaction in acetonitrile-d 3 via in situ 1 H NMR spectrometry; the data are reported in Zenodo (https://zenodo.org/record/3252790). To determine the reaction rate constants for Schemes 1-3 in acetonitrile-d 3 , we use molar concentration data to fit model 1 + 2 + 3 (see ESI †), as shown in Table 1. The estimated parameters at 298 K and their 95% confidence intervals (CI) are presented in Table 3. We specify in Table 4 the NMR peaks used for the fitting and the species associated with each peak in the model. The quality of the model is also summarized in the same table.
At 298 K, the quality of fit for model 1 + 2 + 3 is reasonable, as shown in Fig. 2 and Table 4, and the corresponding PRCR value is in good agreement with the literature value 26 of 94 : 6 and our experimental value of 97 : 3. According to the RMSE and MAPE values, model 1 + 2 + 3 fits the experimental data for each monitored peak well, with slightly elevated values for the naphthol peak (see Table 4), possibly due to the small number of experimental data collected. When estimating the equilibrium constant for Scheme 1, the resulting value is very close to unity (ca. 0.97). Based on this finding and on the pK a and NMR shift values of β-naphthol and the C-alkylated product discussed above, we henceforth set the value of the equilibrium constant equal to 1. By doing so, we significantly decrease the uncertainty in the model parameters while fitting to the experimental data, thus obtaining estimated parameters with significantly narrower confidence intervals. Model 1 + 2 + 3 constitutes a complete representation of the kinetics in acetonitrile. In comparison, the minimal model performs adequately (overall MAPE value of 19%), with larger discrepancies observed for fitting the peaks of naphthol and the C-alkylated product. Based on the kinetic analysis alone, the minimal model may Table 1 Developed models investigated in this study for their suitability to describe the mechanism of the reaction between sodium β-naphthoxide and benzyl bromide in acetonitrile-d 3 and methanol-d 4
Scheme 5 The reaction of benzyl bromide and methanol-d 4 resulting in the formation of benzyl methyl ether.

Kinetics and reaction network in methanol-d 4
With the knowledge obtained from studying the reaction network in acetonitrile, we proceed to study the reaction in methanol-d 4 . We monitor the reaction kinetics in methanold 4 via in situ 1 H NMR spectrometry; the data are reported in Zenodo (https://zenodo.org/record/3252790). The reaction proceeds much more slowly than in acetonitrile, so that it is challenging to estimate the kinetic parameters with precision; on the other hand, this makes data acquisition easier.
In protic solvents like methanol, the selectivity of the reaction towards the O-alkylated product is reduced. 23 As a result of the greater concentration of C-alkylated product, the proton exchange reaction can be expected to have a larger impact on the relative reaction rates and product distribution. The effect on reaction progress can only be quantified by the kinetic modelling. As we will show, the reaction network is more complex in methanol than in acetonitrile and kinetic modelling provides insights into the interplay between different reactions.
We begin by fitting model 1 + 2 + 3 ( Fig. 3(a) and Table 6), which takes into account that sodium β-naphthoxide remains partially sequestered, and thus non-reactive, in the form of β-naphthol (7). The model results in particularly good agreement between the calculated and measured concentrations of total naphthol (1 + 7) and an overall MAPE value of 4.14% is obtained. However, despite the good performance of the model for several peaks, there is a large discrepancy between the calculated and measured amounts of benzyl bromide consumed, suggesting that there is an additional source of consumption for this reactant. This is confirmed by the 1 H NMR spectra in methanol-d 4 where we observe an unexplained peak adjacent to the peak of the C-alkylated product. The existence of a reaction between benzyl bromide and methanol yielding benzyl methyl ether has previously been reported. 42 We confirm this finding by analysing a sample of benzyl bromide in methanol-d 4 , first through a 1 H-13 C HSQC experiment, where we identify the aliphatic CH 2 peak of benzyl methyl ether (4.44 ppm in 1 H NMR, 74.18 ppm in 13 C NMR) in agreement with earlier work; 43,44 and second, through a 1 H-13 C HMBC experiment, where we identify a cross-peak (4.44 ppm in 1 H NMR, 56.09 ppm in 13 C NMR) with consistent carbon chemical shift for the CD 3 group of benzyl methyl ether, illustrating long-range correlations between the proton of the aliphatic CH 2 group with the carbon from the CD 3 group. As a result of these experiments, we include the side reaction between benzyl bromide and methanol in the postulated reaction network in methanol (see Scheme 5) and update our kinetic modelling accordingly (model 1 + 2 + 3 + 5, see ESI †). Experimental concentration data for benzyl methyl ether are derived from the NMR spectra when feasible; peak overlap often occurs due to the proximity of the aliphatic CH 2 peaks monitored for benzyl methyl ether and the C-alkylated product. In some spectra the two peaks are entirely merged, thus we only integrate the spectra in which the two peaks are distinguishable. This results in a different number of fitting data for model 1 + 2 + 3 + 5 than the number of data used for model 1 + 2 + 3 and different concentrations for the C-alkylated product. The parameters of model 1 + 2 + 3 + 5 are regressed from experimental data in methanol-d 4 as listed in Table 5, with concentration profiles shown in Fig. 3(b). Compared to model 1 + 2 + 3, model 1 + 2 + 3 + 5 provides a more complete description of the reaction network and the reduced RMSE and MAPE values for species 1 + 7 and 2 reported in Table 6 indicate that the rates of consumption for benzyl bromide and sodium β-naphthoxide are captured more accurately.
While one might be satisfied with model 1 + 2 + 3 + 5 on the basis of the low errors obtained, analysis of the experimental evidence provides motivation to expand the reaction network further. In their study, Kornblum et al. 23 mention the formation of a third by-product, 1,1-dibenzyl-2-(1H)-naphthalenone (9), a double C-alkylated product resulting from further alkylation of the C-alkylated product with benzyl bromide, as shown in Scheme 4. They report yields of the C-alkylated product and of the double C-alkylated product of 24% and of 4% respectively in ethanol, of 75% and of 10% respectively in 2,2,2-trifluoroethanol, and of 83% and 1% respectively in water. However, in aprotic solvents such as DMF and DMSO, no double C-alkylated The lines correspond to fitting with model 1 + 2 + 3. Symbols denote experimental data. Blue circles correspond to species 1 + 7, red triangles to species 2, green diamonds to species 3 and pink squares to species 4 + 5.

View Article Online
product was observed. 23 In the NMR spectra from our experiments in acetonitrile-d 3 , no peaks corresponding to the double C-alkylated product were observed; thus, the possibility of formation of the double C-alkylated product in acetonitrile-d 3 was dismissed. Given that in ethanol the amount of the double C-alkylated product obtained is appreciable, 23 this product should also be formed in methanol. We thus conduct additional experiments, starting from the deprotonated form of the C-alkylated product reacting with benzyl bromide, to confirm the formation of the double C-alkylated product, as shown in Scheme 4. We also confirm that the protonated form of the C-alkylated product is incapable of further alkylation, as no sign of the double C-alkylated product is observed in 1 H NMR spectra obtained after leaving the C-alkylated product with benzyl bromide in a vial for 24 hours (in the absence of any deprotonating agent, e.g., a base). This finding is consistent with the fact that naphthol, the protonated form of β-naphthoxide, is also incapable of reacting with benzyl bromide; the two species (naphthol and the C-alkylated product) have a very similar chemical structure, and thus similar reactivity. On the basis of these findings, we add Scheme 4 to the reaction network.
Model 1 + 2 + 3 + 4 + 5 incorporates Scheme 4 and comprises all the knowledge obtained from experimental and kinetic modelling observations. Concentration profiles for model 1 + 2 + 3 + 4 + 5 at 298 K are shown in Fig. 3(c), with the estimated parameter values listed in Table 5. As shown in Table 6, the PRCR value for model 1 + 2 + 3 + 4 + 5 of 73 : 27 is in very good agreement with our experimental value of 72 : 28, offering a significant improvement over the values found with model 1 + 2 + 3. A similarly good performance is View Article Online achieved with model 1 + 2 + 3 + 5, with a PRCR of 75 : 25, due to the slow rate of the reaction producing the double C-alkylated product. However, this model can be expected to deviate further from the experimental data as time progresses and the additional loss of C-alkylated product occurs. The same RMSE value is obtained for model 1 + 2 + 3 + 5 and model 1 + 2 + 3 + 4 + 5 (0.0008 mol dm −3 ) when the same data for the four main peaks are considered ("Restricted average" column in Table 6). RMSE values for each model when all data are considered ("Overall" column in the Table 6) confirm the superiority of model 1 + 2 + 3 + 4 + 5.

Comparison of kinetics in both solvents
Overall, based on the reaction-rate constants listed in Tables 3 and 5, O-alkylation proceeds faster than C-alkylation in both solvents, as indicated by the larger k O values and the resulting higher selectivities to the O-alkylated product. By comparing the values shown in Tables 3 and 5, we observe that the O-alkylation pathway proceeds much faster in the polar aprotic solvent than in the polar protic solvent. At 298 K, the estimated value of k O is nearly two orders of magnitude larger in acetonitrile than in methanol. In both solvents, the reaction-rate constants for C-alkylation, k C , are of a similar order of magnitude. The reaction constant for double C-alkylation, k DC , is found to be very similar to that of C-alkylation in methanol. The formation of benzyl methyl ether proceeds at a much slower rate than other reactions, so this only has a small impact on the observed selectivity. The measured PRCR of 97 : 3 observed in acetonitrile demonstrates a clear preference towards the O-alkylated product, in agreement with previous studies in polar aprotic solvents; 23,25,26,28 in particular, the value of 97 : 3 at 298 K obtained here is in excellent agreement with the value of 94 : Table 5 Estimated reaction-rate constants k O (dm 3 mol −1 s −1 ) and k C (dm 3 mol −1 s −1 ) for Scheme 1, equilibrium constant K eq (−) for Scheme 3, reaction-rate constant k BME ·[8] (s −1 ) for Scheme 5, and reaction-rate constant k DC (dm 3 mol −1 s −1 ) for Scheme 4, with their 95% confidence intervals (CI), obtained based on the various models developed for methanol-d 4   a This is the calculated PRCR given by eqn (4). b This is the calculated PRCR given by eqn (6)

Temperature dependence of the reaction rate constant
We conduct further analysis using the Arrhenius equation to estimate the pre-exponential factor (A) and the activation energy (E a ), and the Eyring equation to obtain the enthalpy of activation (Δ ‡ H _°) and entropy of activation (Δ ‡ S _°). We use the experimental concentration profile data obtained at three temperatures in acetonitrile (at 298, 313 and 323 K) and in methanol (at 298, 313 K and 328 K) in the parameter estimation. This corresponds to 136 data points in acetonitrile and 295 data points in methanol. Due to the narrow temperature range studied, the thermodynamic quantities are assumed to be temperature-independent. Model 1 + 2 + 3 is used for acetonitrile and model 1 + 2 + 3 + 4 + 5 for methanol. The estimated thermodynamic parameters obtained are listed in Tables 7 and 8 for acetonitrile and methanol, respectively. The quality of fit across all three temperatures and both solvents is found to be highly satisfactory: the Arrhenius model gives an RMSE of 0.0030 mol dm −3 in acetonitrile and 0.0027 mol dm −3 in methanol, while the Eyring model yields RMSE values of 0.0022 mol dm −3 in acetonitrile and 0.0027 mol dm −3 in methanol. A more detailed analysis for each peak can be found in the ESI † (Fig. S1-S12 and Tables S6-S10).
Of particular interest is the standard-state Gibbs freeenergy difference between the C-and the O-alkylation transition-state structures in a given solvent, at a specific temperature (ΔΔ ‡ G̲°; L C−O T ð Þ), which can be computed as the difference in the standard-state free energies of activation of the C-and O-alkylation pathways (cf., eqn (13)). This quantity dictates the selectivity of the reaction at a specific temperature and is related to the reaction-rate constant ratio, which can be obtained from eqn (12) assuming the transmission coefficient to be unity. Based on the definition of ΔΔ ‡ G̲°; L C−O T ð Þ, a positive sign indicates a selectivity preference towards the O-alkylation pathway at a specific temperature, whereas a negative sign indicates a selectivity preference towards the C-alkylated product.
Thermodynamic analysis based on our experimental data results in ΔΔ ‡ G̲°; L C−O 298 K ð Þ values of 7.0 kJ mol −1 in acetonitrile and of 2.8 kJ mol −1 in methanol, corresponding to PRCR values for O-alkylated product to C-alkylated product of 97 : 3 and 72 : 28, respectively. This provides a quantification of the major impact of the solvent on the relative reaction barriers. To gain further understanding of the way in which the solvent affects the reaction, we turn to electronic-structure methods and assess the relative stabilisation of the reaction species.

Reaction analysis with quantum mechanical calculations
Density functional theory (DFT) methods have been previously employed for predicting reaction kinetics, 2,45-47 with promising outcomes. In our study, we use computational chemistry methods to investigate the mechanism and kinetics of Scheme 1 for the first time, focusing particularly on the transition states and reaction energetics for Scheme 1.
3.7.1 Transition state structures. We locate the transition state structures of the two alkylation pathways in vacuum following the hypothesis of Kornblum et al. 23 that the O-alkylation pathway proceeds via a transition state with a linear arrangement of the oxygen, carbon and bromine atoms, while the C-alkylation pathway proceeds via a linear arrangement between the nucleophilic carbon, the electrophilic carbon and the bromine atom. The transitionstate structures obtained for the O-and C-alkylation pathway in vacuum are shown in Fig. 4. Structural parameters and frequency analysis for the optimized transition-state structures are reported in Tables S11-S19 in the ESI. † For the O-alkylation pathway, two transition states are identified, TS O1 and TS O2 , that differ in the orientation of the plane of the ring of the β-naphthoxide ion with respect to the plane of the ring of benzyl bromide. In TS O1 the planes of the two aromatic ring systems are perpendicular, with a dihedral angle of 88°, while in TS O2 the two planes are parallel, with a dihedral angle of 178°. TS O1 resembles closely the structure postulated by Kornblum et al.; 23 however, TS O2 has a lower electronic energy by 5 kJ mol −1 (all levels of theory result in very similar geometries; the structures shown in Fig. 4 were optimized at B3LYP/6-31+G(d) followed by a single point energy calculation at MP2/cc-pVDZ). Thus, this pathway makes a much greater contribution to the reaction rate. We find that TS O2 is always of lower energy than TS O1 when the structures are optimized in acetonitrile and methanol. As a result, we assume hereafter that the O-alkylation pathway proceeds only via TS O2 . In the C-alkylation transition-state structure, in line with what was suggested by Kornblum et al., 23 the carbon atom (C 3 ) from which the bromine atom departs is above (or below) the plane of the β-naphthoxide ion, which attacks the carbon with the departing bromine atom rearward, as is consistent with an S N 2 mechanism. The C-alkylation pathway thus proceeds via a quasi six-atom transition state. The main difference between the structure we find here and the one suggested by Kornblum et al. 23 is that the sodium cation stabilizes the benzene ring instead of being close to the bromide ion. The stabilization occurs through cation-π interaction (non-covalent), keeping the sodium cation away from the centre of the benzene ring at a distance of 2.54 Å in vacuum and of 2.55 Å in methanol, and at a slightly longer distance in acetonitrile (2.75 Å). According to Ma and Dougherty 48 and Amicangelo and Armentrout, 49 the free energy of binding between benzene and sodium cation is 27 kcal mol −1 in the gas phase and decreases as we move to larger cations, e.g., K + , NH 4 + , etc. This trend supports the idea of coulombic forces holding the key role in the specific interaction's strength. Our study provides structural details on the transition-state structures of the two pathways that would be impossible to discern from experimental studies alone. We also gain valuable insights into mechanistic details, for example the orientation of the naphthalene ring and the involvement of the cation in the transition-state structures of the two pathways, that complement our knowledge of the mechanism of Scheme 1. 3.7.2 Relative energy barriers. Next, we investigate the extent to which DFT calculations can predict relative energetic quantities, focusing initially on the relative activation barrier between the two pathways in Scheme 1, values. Among all levels of theory tested, M05-2X/6-31G(d) yields the best agreement with experiments in both solvents, predicting a ΔΔ ‡ G̲°; L C−O 298 ð Þ value of 6.67 kJ mol −1 in acetonitrile (expt. 7.0 kJ mol −1 ) and 1.93 kJ mol −1 in methanol (expt. 2.8 kJ mol −1 ). The free-energy of solvation difference between the C-and the O-alkylation transition-state structures is one of the two key factors for accurately predicting the value of ΔΔ ‡ G̲°; L C−O T ð Þ (the other being the gasphase free-energy difference). The good performance of M05-2X/6-31G(d) in predicting ΔΔ ‡ G̲°; L C−O 298 ð Þ could perhaps be anticipated from its mean unsigned error (MUE) of 2.68 kJ mol −1 for non-aqueous solvation free energies (where methanol is excluded) on a test set of 2072 neutral solutes. 36 Indeed, several authors have discussed the need for adhering to the computational protocol used in the parameterization of each continuum solvation model as a primary requirement in achieving quantitative results using this model. 45,50,51 The accuracy of M05-2X/6-31G(d) was the highest among all the levels of theory used.
Specifically for acetonitrile, Marenich et al. 36 reported that solvation free-energy calculations at the M05-2X/6-31G(d) level of theory result in a MUE of 2.64 kJ mol −1 , based on a sample of 7 neutral solutes. The reported error value increased to 24.70 kJ mol −1 when a sample of 69 ionic solutes was considered. Despite the large size of the sample of ionic solutes, the error value is appreciable and gives an indication of the level of accuracy expected in the solvation free-energy calculations with the SMD model when charged species are We further report QM calculations for the solvation free energy (ΔG _°, solv ) for the reactants and the two transitionstate structures using the 4 levels of theory in common between this study and the study of Marenich et al.; 36 they can be found in Table S22 and Fig. S16 in the ESI. † We focus our analysis on the M05-2X/6-31G(d) level of theory, given its good performance in calculating relative free energies. A comparison of the solvation free-energy values in the two solvents indicates that benzyl bromide and the transition-state structure for the O-alkylation pathway are destabilized in moving from acetonitrile to methanol, exhibiting a lower (more negative) free-energy of solvation in acetonitrile than in methanol, with absolute energy differences of 3.56 kJ mol −1 for benzyl bromide and 4.41 kJ mol −1 for the transition state. On the other hand, sodium β-naphthoxide and, to a smaller extent, the transition-state structure for the C-alkylation pathway are subject to a decrease in solvation free-energy from acetonitrile to methanol (i.e., stabilisation) with absolute energy differences of 4.97 kJ mol −1 and 0.33 kJ mol −1 , respectively. Overall, the two reactants undergo a net stabilisation via a decrease of 1.41 kJ mol −1 in solvation free energy while the two transition states undergo a destabilisation via an increase of 4.08 kJ mol −1 in solvation free energy. The main driver for the change in selectivity from acetonitrile to methanol is due to the destabilisation of the transition-state for the O-alkylation pathway relative to the C-alkylation transitionstate, which is barely affected.
3.7.3 Computed reaction rate constants. We also compute absolute reaction rate constants for the O-and the C-alkylation pathways using eqn (9) for each level of theory tested (listed in Tables S20 and S21 in ESI †). The limitations of the G _ CDS term in capturing the short-range nonelectrostatic interactions may be compensated by the increasingly accurate prediction of the bulk electrostatic interactions, which is again subject to the formalism and parameterization of the functional used. 52 As a result, different levels of theory capture better either relative energy differences or absolute energies. For example, in our study, M05-2X/6-31G(d) has the best agreement with experiments when computing ΔΔ ‡ G̲°; L C−O 298 ð Þ values (MAPE = 18%) but fails to capture absolute reaction rate-constant values by two orders of magnitude. There are, of course, levels of theory, such as B3LYP/6-31+G(d) and M05-2X/cc-pVTZ, that perform reasonably well for both relative and absolute reaction metrics. But generally, when computing absolute reaction rate constants for the O-and the C-alkylation pathway, we resort to functionals with more extensive basis sets, such as M05-2X/cc-pVTZ, which will provide us with more accurate estimations regarding bulk electrostatic interactions. In our study, the best agreement with experimental reaction rate constants (on average) for the two alkylation pathways is achieved in acetonitrile by M05-2X/cc-pVTZ (MAPE = 63%) and in methanol by B3LYP/6-31+G(d) (MAPE = 29%). As might be expected, augmenting the 6-31G(d) basis set with a diffuse function (+) to account for the presence of the bromide anion improves the functional's performance. As such, B3LYP/6-31+G(d) provides the best overall agreement with measured reaction rate constants for both solvents (MAPE = 47%), with a correct prediction for selectivity preference of the reaction (MAPE = 47%), also in both solvents.
We use the most accurate QM-computed rate constants instead of the fitted rate constants to simulate the reaction in both solvents and to assess the impact of the errors in the QM-computed rate constants. For the side reactions, we use the experimentally-fitted reaction constants. When this is done using experimental reaction conditions, the 63% error in the kinetic constants in acetonitrile translates to a RMSE value of 0.0120 mol dm −3 and a MAPE value of 81.5%. In methanol, where there are more side reactions, the 29% error in the rate constant translates into an RMSE value of 0.0043 mol dm −1 , a restricted average RMSE value of 0.0051 mol dm −3 , a MAPE value of 29% and a restricted average MAPE of 17%. These findings indicate the large extent to which errors in rate constants propagate as errors in the concentrations. This, combined with the variability of the accuracy of each level of theory with solvent, highlights the need for further development of solvation models.

Complementarity of investigative methods
It is useful to reflect upon the roles played by experiments, kinetic modelling and DFT calculations in elucidating and quantifying the mechanism and kinetics of the Williamson reaction studied here. Clearly, experimental studies via quantitative 1 H NMR, 1 H-13 C HSQC, 1 H-13 C HMBC, and 13 C NMR all play a vital role in understanding which species are present and the evolution of the concentrations of several species as a function of time, in different solvents and at different temperatures. This is vital in guiding modelling efforts, be it kinetic or DFT modelling.

View Article Online
The kinetic modelling plays several roles. On the one hand, it guides the interpretation of experimental data, forcing its re-evaluation when there is a mismatch between data and model predictions. For example, the deterioration in the quality of fit in benzyl bromide concentration observed from model 1 to model 1 + 3 (addition of the proton exchange reaction) in methanol helps to focus the analysis on a different reaction involving benzyl bromide, despite the small size of the peaks for benzyl methyl ether. On the other hand, the derivation of the reaction rate constant values provides a quantitative understanding of what matters for the reaction network. The proton exchange reaction is found to have similar equilibrium constants in both solvents, but only has a small impact on reaction progress in acetonitrile. The impact is much more significant in methanol, because of the smaller reaction rate of O-alkylation and the higher (detectable) reaction rate of double C-alkylated product formation. Finally, the kinetic models developed make it possible to predict reaction progress for different starting concentrations and at different temperatures than those used in the experimental studies. It is essential to capture in this way our understanding of the reaction networks so that this can be used in the design of further experiments or in reactor design calculations.
DFT calculations enable us to achieve a more detailed understanding of the reaction mechanism. Computations of structures appear reliable and help to test existing postulates on the geometries of the transition states, providing insights into the differences between the two pathways. Our experience here shows that those levels of theory that were used to parameterize the solvation model provide a good level of accuracy in terms of relative energy values, which are mostly found to be well within chemical accuracy (1 kcal mol −1 or 4.18 kJ mol −1 ), with the exception for M05-2X/6-31+G(d) in methanol. Such accuracy is sufficient to ascertain that the O-alkylation pathway is always more favoured over the C-alkylation pathway in acetonitrile than in methanol. This is not enough to achieve quantitative accuracy when we calculate the rate constants and relative rates; hence, the kinetic modelling and experimental data are useful in identifying the most reliable levels of theory for this reaction. The analysis nevertheless indicates that switching solvents from acetonitrile to methanol leads to the destabilisation of the transition state structure for the O-alkylation pathway and a consequent change in regioselectivity as the C-alkylation reaction becomes more competitive.

Conclusions
In this work, we have pursued a three-faceted approach that integrates experimental probing, kinetic modelling and quantum-mechanical calculations to provide insights into the mechanism of a Williamson ether synthesis reaction. Each methodology contributes towards the elucidation of mechanistic aspects as well as providing validation of the observations made with other approaches. Specifically for the regioselective Williamson reaction between sodium β-naphthoxide and benzyl bromide investigated here, the findings from a three-faceted analysis for solvents acetonitrile and methanol have brought to light subtle changes in product distributions and the balance between reaction pathways. Detailed experimental work in the two solvents has revealed an underlying equilibrium between two of the reaction species, one of the reactants and one of the products, taking place in both solvents, as well as two additional reactions occurring in the presence of methanol: the further reaction of the C-alkylation product to form an additional third product, and a side reaction occurring between methanol and one of the reactants.
The complexity brought about by the change of solvent medium highlights the role of kinetic modelling as an integral component of the mechanistic analysis. Different kinetic models were developed to account for the species observed experimentally. The qualitative and quantitative description of the experimental concentration profiles of the reaction species as well as the statistical significance of the estimated parameters were carefully evaluated. Discrepancies between experimental data and kinetic model led to the reevaluation of experiments and incorporation of the newly obtained knowledge in the kinetic models. Model 1 + 2 + 3 for acetonitrile and model 1 + 2 + 3 + 4 + 5 for methanol were found to give the best description of reaction network and kinetics with RMSE values of 0.0028 mol dm −3 and 0.0008 mol dm −3 (for the main 4 reaction components), respectively. Estimated thermodynamic parameters were also obtained based on the corresponding model for each solvent. For the quantity ΔΔ ‡ G̲°; L C−O T ð Þ that dictates the reaction selectivity at 298 K, values of 7.0 kJ mol −1 in acetonitrile and of 2.8 kJ mol −1 in methanol were obtained, corresponding to selectivity ratio values for O-alkylated over C-alkylated product of 97 : 3 and 72 : 28, respectively. Quantummechanical calculations provided mechanistic insights into the two reaction pathways such as the identification of transition-state structures for the O-and C-alkylation pathways, as well as an alternative transition-state structure for the O-alkylation pathway as opposed to the existing literature hypothesis. Further understanding of the way in which the solvent affects the reaction was obtained by assessing the relative stabilisation of the reaction species in the two solvents through calculation of the liquid free-energy and free-energy of solvation, using levels of theory consistent with the parametrization of the solvation model SMD. The solvent does not appear to alter the transition states for the two pathways, but rather acts via stabilization/destabilization of the reaction species: the O-alkylation transition-state structure and benzyl bromide are always more stable in acetonitrile than in methanol (maximum stabilization difference observed for G _°, L of 6.3 kJ mol −1 at the M05-2X/6-31G(d) level of theory and 3.8 kJ mol −1 at the M05-2X/6-31G +(d,p) level of theory), whereas the C-alkylation transitionstate structure and sodium β-naphthoxide are always more stable in methanol than in acetonitrile (maximum View Article Online stabilization difference observed for G _°, L of 1.2 kJ mol −1 at the M05-2X/6-31+G(d,p) level of theory and for ΔG _°, solv of 5.5 kJ mol −1 at the same level of theory). Within the reported accuracy of the computational methods used, quantitative conclusions were also drawn that indicate that the O-alkylation pathway was always favoured in acetonitrile.
No definite conclusions could be drawn in methanol due to smaller energy difference between the two pathways. The error trend was the same in both solvents, showing that the degree of overestimation or underestimation is consistent for each computational model. A range of values can be obtained for estimating Δ ‡ G̲°; L C−O 298 ð Þ with a MUE of 3.52 kJ mol −1 in acetonitrile and 4.02 kJ mol −1 in methanol.
The findings in our work can be used to inform investigations of solvent effects on Williamson ether syntheses and other organic reactions. For the specific scheme considered here, it is possible to use the most accurate levels of theory identified to conduct a preliminary computational assessment of selectivity and rate in a range of other solvents, thereby guiding any experimental programme. For other substrates, similar trends can be expected for the effect of solvents on kinetics. The proposed transition state structures can be adapted for new substrates and used as a starting point to accelerate the search for a detailed mechanistic understanding. More broadly, our work provides a blueprint for the investigation of solvent effects in organic reactions, highlighting the value and complementarity of different analysis tools in such studies and demonstrating specific types of solvent effect on mechanism and kinetics (e.g., proton exchange, structural changes in transition state, reaction between reactant and solvent). Our systematic analysis of the accuracy of levels of theory shows that QM calculations can play an important role in the qualitative understanding of mechanistic impacts, but that the choice of an appropriate level of theory for quantitative analysis must be validated with experimental data.
Overall, this work highlights the value of a three-faceted approach to elucidate the impact of reaction conditions on the reaction outcome. Recognising that measurements, kinetic modelling and quantum mechanical calculations are associated with different types of uncertainty, the combination of these three methods of investigation offers a robust way to examine complex phenomena such as regioselectivity: in some cases, findings can be corroborated by two methods, providing added confidence to draw conclusions, while other aspects (e.g., the geometry of transition states) can only be probed with one of the methods. In practice, the three-faceted approach encourages the iterative improvement of the understanding of the system under study as initial hypotheses are challenged by using different types of analysis and must be revised. The principles of the three-faceted approach are generic, so that it can be applied to other reaction studies where quantitative analysis and prediction of quantitative chemical concepts are needed. The use of the proposed approach is particularly promising for the analysis of reactions that involve expensive or scarce components or that can result in toxic side products. The approach could be applied to a small number of solvents to identify a QM strategy (e.g., Struebing et al. 45 ) to screen other solvents computationally, reducing the number of experiments required.

Data statement
Data underlying this article may be accessed on Zenodo at 10.5281/zenodo.3252790 and used under the Creative Commons Attribution licence.

Conflicts of interest
There are no conflicts of interest to declare.