Neutron Compton Scattering: from proton momentum distribution to muonium hyperfine coupling constant in the isopropyl radical

We establish a fast and reliable benchmarking protocol for predictions of Muon Spin Resonance observables. To this end, we apply neutron Compton scattering (NCS) to study the nuclear momentum distributions of the proton and deuteron in the condensed phase of the isopropyl and d-isopropyl alcohols. By subtracting the time-of-flight NCS spectra of both compounds we demonstrate that the proton momentum distribution in the OH group of isopropanol and the deuteron momentum distribution in the OD group in d-isopropyl can be studied selectively. The site-selective application of the NCS method enables the calculation of the magnitude of the frequency isotope effect for the proton in OH along the hydrogen bond direction. By comparing the magnitude of the frequency isotope effect with values predicted for simple model potentials we are able to perform the appraisal of the degree of anharmonicity of the OH proton environment. Assuming that the effective potential felt by the OH proton along the hydrogen-bond direction can be satisfactorily described by the Morse potential, we are able to calculate its dissociation constant D and decay constant a. Finally, assuming that the same Morse potential describes the local binding of Muonium in the mioniated isopropyl radical, we are able to predict its width of momentum and position distributions and the kinetic and zero-point energy. Based on these results, we are able to provide a conservative bound for the magnitude of the isotope effect on the muonium hyperfine interaction without resorting to a complicated and computationally expensive methodology based on the application of path integrals.


Muon spin rotation and its challenges
Materials science at the nanoscale over last few decades has profited from (and largely owed its success to) the development of a few very special and unique spectroscopic techniques that do not have the stigma of being volume-averaged. One of the most prominent examples of such a development owing to the site-specific method sensitivity is the muon spin resonance (μSR) [1]. The supremacy of μSR over volume-sensitive nuclear magnetic resonance (NMR) is best manifested in the spin spectroscopy of the solid-state where the interpretation of μSR spectra is greatly simplified by the absence of spectral line broadening due to the lack of the long-range spin-spin couplings (otherwise present in the NMR spectra) thanks to the local detection character of the muon, acting as a local magnetic moment probe [2]. Namely, instead of a static NMR coil wound around a macroscopic volume of a sample, one has a fast-moving 'coil-projectile', the muon particle, picking up signals from local magnetic environments. During the muon's 2.2 μs lifetime, muonium can enter into compounds such as muonium chloride (MuCl) or sodium muonide (NaMu) [3]. Due to the mass difference between the antimuon and the electron, muonium (μ+e−, hereinafter Mu) is more similar to atomic hydrogen (p+e−) than positronium (e+e−) [3]. Its Bohr radius and ionization energy are within 0.5% of hydrogen, deuteron, and tritium, and thus it can usefully be considered as an exotic light isotope of hydrogen with a mass of l/9 of H [3].
The modelling of Mu in condensed matter systems up to date has dealt with properties such as Mu formation in gases, liquids, and solid insulators, and the rate constants of its reactions with molecules [2]. Despite the unquestioned success in the experimental department, also μSR is not entirely free from problems as far as ab initio modelling is concerned [3]. The main challenge here is due to the fact that the Zeeman interaction of the unpaired electron, the muon and of additional nuclei, as well as the electron-muon and electron-nuclear hyperfine interaction, lead to very complicated μSR spectra [4]. The variation of the muon-electron coupling constants with radical structures is similar to that of the analogous known proton-electron coupling constants. However, the small muon mass greatly amplifies any vibrational isotopic effect in these coupling constants as well as all other properties dependent on the vibrations of the radical. One could anticipate that for muonium adduct molecules, the isotope effect is so large that the Born-Oppenheimer approximation (BOA) is invalid, a situation certainly rarely encountered in the case of ab initio modelling of electronic and nuclear degrees of freedom in condensed-matter systems. For isotopomers of the H 2 + molecule Webster and McKenna showed by a series of non-adiabatic calculations beyond the Born-Oppenheimer approximation, that the key isotope effect arises through the mass dependence of the zero-point vibrational energy [5]. On the whole, one has to anticipate that the Mu-adduct species sample points on a property surface well distant from the equilibrium position. For example, a variation-perturbation approach, applied to investigate isotope-dependent zero-point vibrational corrections to one-electron properties of a MuMu adduct water molecule in its ground electronic state, has shown that the root-mean-square displacement involving the O-Mu internuclear distance is ca. 11 pm with a root-mean-square displacement in the internuclear angle for MuOH of 12.73 degrees and MuO of 14.71 degrees [5]. Ab initio modelling under such circumstances is very difficult, rendering most DFT methods treating nuclei classically within the BOA completely unusable. This can be demonstrated owing to the recent rapid development of Path-Integral Molecular Dynamics (PIMD) modelling methods [6]. For example, the on-the-fly ab initio density functional PIMD simulations, which can account for both the nuclear quantum effect and thermal effect, when applied to the evaluation of the structures and isotropic hyperfine coupling constants (HFCCs) for muoniated and hydrogenated acetone radicals in vacuo, exposes dramatic deficiencies of the DFT methodology with classical nuclei. The HFCC value from a simple geometry optimization calculation without both the nuclear quantum effect and thermal effect is 8. 18 MHz, and that by standard ab initio molecular dynamics simulation with only the thermal effect and without the nuclear quantum effect is 0.33 MHz at 300 K, where these two methods cannot distinguish the difference between muoniated and hydrogenated acetone radicals. In contrast, the reduced HFCC value of the muoniated acetone radical by the PIMD simulation is 32.1 MHz, which is about 8 times larger than that for the hydrogenated radical of 3.97 MHz with the same level of calculation [6].

The need for a mass and site-selective dynamic probe
The calculation of many properties of muoniated radicals starts from the appraisal of the extent of quantum delocalisation of Mu, both in momentum and position representation. In general, the momentum distribution can be obtained from the Fourier transform of the end-to-end-distribution of the linear (or 'open') polymer PIMD calculation, but not from the ring (or 'closed') polymer PIMD. Moreover, the PIMD methods using generalised Langevin equation (PI-GLE) approach employ a 'trick' for avoiding the need for a large number of replicas, as far as the physical temperature of the system is low. However, it is only by using this 'artificial noise' thermostat with adequate and well-tested parameters that one can create random fluctuations in the motion of the nuclei which are similar to the zero-point broadening of the real quantum delocalization. To speed up the calculations and base them on the closed-path PIMD scheme, one might resort to an additional approximation that the nuclear wave function is node-free so that it can be obtained as the square root of the nuclear density in coordinate-space, or the square root of the Fourier transform of nuclear momentum distribution. However, despite the sophistication of the PIMD-based methodology, the results calculated in vacuo will always be far off the real values measured in bulk systems as the explicit presence of solvent molecules has a major effect on in vacuo ab initio calculations for the quantitative reproduction of experimental data. Clearly, there is a need for a quick and model-free method for the appraisal of properties of muoniated radicals under realistic bulk conditions. In this context, the employment of the Deep Inelastic (DINS), or neutron Compton Scattering (NCS) seems a plausible choice. NCS, as the only neutron-based spectroscopic technique, offers true mass selectivity by recording a combination of Doppler-broadened recoil peaks for different nuclei present in the sample of interest [7][8][9][10][11][12][13][14][15][16][17][18][19]. Moreover, in the limit of the incoherent and impulse approximations, valid in NCS, the width of a Doppler-broadened NCS peak of a given nucleus is proportional to its nuclear momentum distribution, n(p), being the square of the absolute value of the Fourier transform of the spatial wave function [7][8][9][10][11][12][13][14][15][16][17][18][19]. Importantly, the small mass of the proton/deuteron (H/D) leads to an NCS recoil peak that is well separated from the peaks of heavier nuclei, which appear as nearly elastic contributions. This feature, in combination with the fact that NCS is sensitive to total bound scattering cross-sections, makes muoniated radicals containing H and D ideal candidates for the NCS measurements.
This contribution seeks to pave the way towards establishing a robust model-free protocol for the appraisal of basic properties of the muoniated radicals and a reliable benchmarking protocol for the PIMD-based predictions. As a feasibility study of the newly proposed computational protocol, we analyse here the results of NCS measurements on two model bulk systems with a similar local environment around the OH group, the isopropyl and d-isopropyl alcohol (figure 1).
The idea to perform a tandem NCS measurement of a native and deuterated counterpart of the same molecular system, both differing only in the presence (or absence) of H(D) in the OH group, renders an otherwise mass-selective technique also site-selective. Namely, by virtue of the mass-selectivity, the direct NCS measurement of the deuterated compound yields the characterisation of the deuteron momentum distribution in the OD group. Furthermore, by subtracting the NCS spectra of both compounds one obtains the selective spectrum of the proton momentum distribution in the OH group of the isopropanol molecule that can then be compared with the deuteron momentum distribution in the OD group in d-isopropyl. This site-selective application of the NCS method enables the calculation of the magnitude of the frequency isotope effect for the proton in OH along the hydrogen bond direction. By comparing the magnitude of the frequency isotope effect with values predicted for simple model potentials we are able to perform the appraisal of the degree of anharmonicity of the OH proton environment. Assuming that the effective potential felt by the OH proton along the hydrogen bond direction can be satisfactorily described by the Morse potential, we are able to calculate its dissociation constant D and decay constant a. Finally, assuming that the same Morse potential describes the local binding of Muonium in the mioniated isopropyl radical, we are able to predict its width of momentum and position distributions and the kinetic and zero-point energy. Moreover, based on the latter results we provide a conservative bound for the magnitude of the isotope effect for the hyperfine interaction of the muonium in the muoniated isopropyl in the condensed phase, all without resorting to a complicated and computationally expensive methodology based on the application of path integrals.

Materials and methods
We have performed neutron experiments on isopropanol and isopropanol-2-d 1 (hereafter referred to as d-isopropanol), at 250 K and 300 K, by concurrently applying the mass-selective NCS (MANSE), enabling distinct separation of H and/or D peaks, and neutron transmission (NT) at the VESUVIO+beamline in the Target Station 1 (TS1) of the ISIS pulsed neutron and muon source at the Rutherford Appleton Laboratory, UK [10,14,[20][21][22][23][24][25][26][27][28][29]. Isopropanol (99.0%) and d-isopropanol (98.0%) neat liquid samples were purchased from Sigma-Aldrich. Two square-shape (90×90 mm 2 ) aluminium containers for liquid samples were used for the experiments: one of 0.5 mm thickness for the isopropanol sample, and one of 1.8 mm thickness for the d-isopropanol sample. Both liquids were loaded into the containers through a pair of small inlets in the containers front walls using a set of syringes. Following this, the inlets were sealed. The empty containers were then separately measured at the same temperatures as the samples. Both sample and empty-container neutron transmissions were measured as a function of the incident neutron energies, normalising to the transmission spectrum of the empty instrument, a protocol routinely applied in broadband transmission experiments on VESUVIO [30][31][32][33][34]. The mass of the isopropanol sample was 0.89 g. Given the total bound scattering cross-section values of H, C, and O and the sample stoichiometry, C 3 H 8 O, the mass of 0.89 g corresponded to the sample transmission for epithermal neutrons (above the incident neutron energy of 1 eV) at the level of 0.86. The mass of the d-isopropanol sample was 1.58 g, which, given the sample stoichiometry, C 3 H 7 OD, resulted in the sample transmission for epithermal neutrons at the level of 0.82.
The NCS data treatment has been described in great detail elsewhere [9,10,35] and thus only the main points will be mentioned here. In essence, in the Impulse Approximation [36] limit, a MANSE spectrum measured by a particular detector (at a given scattering angle θ) and at a given time-of-flight (TOF) value, t, is a sum of neutron Compton recoil peaks. Mathematically, [9,10] 7. after the subtraction of the backgrounds, a final fit of peak widths was performed and the final averaged values of the peak widths and relative peak intensities were obtained; 8. fitted curves of the recoil peaks from masses other than the proton and the fitted curves representing the FSE contributions to the proton recoil peaks (see the second term in equation (2)) were subtracted and the corrected TOF NCS data; 9. the isolated TOF NCS proton recoil peaks were focused in the longitudinal momentum space of the proton (hereafter referred to as the proton Y space data) for precision fitting of the proton NMD using purely Gaussian expression (expressed as the first part of equation (2)); 10. in the case of the TOF NCS deuteron recoil peaks, no longitudinal momentum space focusing was performed and the final values of the D peak widths were obtained as averages over sequential, detector-bydetector, fitting of the TOF NCS spectra containing D recoil peaks; 11. in order to obtain the proton momentum distribution of the proton in OH, a purely Gaussian fit of the difference between the focussed proton Y space data of isopropanol and d-isopropanol was performed. Before subtraction, the proton Y space data of d-isopropanol were scaled by the factor of 7/8 reflecting the ratio of the number of protons in the formula units of both systems.
Before we proceed to the detailed description of the main steps of NCS data treatment, outlined in points 1-11 above, illustrated in figure 2, one important remark is in order. As described in points 9-11 of the protocol, the recoil peak of hydrogen and deuteron undergo a completely different treatment. The deuteron recoil peaks are fitted sequentially, detector-by-detector, and the values of the fitted widths are recorded in order to be averaged over the entire range of detectors used in the data analysis. The proton recoil peaks, however, after being provisionally fitted on the detector-by-detector basis, undergo further treatment. They are isolated in the TOF domain, transformed into the longitudinal momentum domain and then focused and rebinned in this latter domain for higher precision fitting. Due to this substantial difference in the treatment, the proton and deuteron recoil peak fitting cannot be easily illustrated in on simple graph. Mathematically, the focusing of proton recoil peaks in the proton longitudinal momentum space without their prior isolation (and thus retaining the deuteron recoil peaks in the focused spectra) would always result in smearing out the deuteron recoil peaks, a kinematic effect due to different loci of the proton and deuteron recoil lines for each scattering angle (detector). Converesely, fitting the TOF NSC spectra sequentially, detector-by-detector, does not allow for the longitudinal momentum space focusing of either of the two recoiling nuclear species. Thus, as an illustration, in figure 2 the TOF NCS data treatment scheme is shown only for one chosen VESUVIO detector. Such choice retains the possibility of graphically showing the logic behind the data treatment with one drawback, however. Namely, the procedure of the focusing of the proton NMD in the proton longitudinal momentum space is illustrated only as a mathematical transformation without showing their added benefit of increasing the signal-to-noise ratio of the proton recoil peak data and thus the precision of the proton NMD fit. Figure 2 shows data recorded for isopropanol and d-isopropanol at the temperature of 300 K. The data were recorded by the detector placed at a distance of 0.52 metres from the sample position in the forward scattering direction at a scattering angle of 63.4 degrees. The position and scattering angle values pertaining to this detector were obtained via a separate calibration procedure described elsewhere [22,23]. The TOF NCS data, prior to further analysis, were rebinned between 110 and 385 microseconds with a bin width of 3 microseconds. The lower binning boundary was set to a relatively high value of 110 microseconds, a procedure routinely applied in NCS data reduction in order to avoid the contamination of the TOF NCS signals with spurious signals from resonant neutron absorption in the Au foils of the VESUVIO detection setup [25][26][27]. It is worth noting that this procedure does not affect the precision of the fitting of the width of the proton NMD. The reason for this is twofold. Firstly, the width of a Gaussian profile can be very robustly estimated from the fits of the recoil peaks as most of the peak profile is anyway retained between 110 and 385 microseconds. Secondly, the fit precision is augmented by the presence of an additional physical constraint, the requirement that an NMD be an even function of the longitudinal momentum distribution variable Thus, in principle, it is enough to record one shoulder of an NMD and then apply symmetrization in order to recover the full profile [9-11, 14, 15, 18, 19, 43-45].
The top and middle panels of figure 2 show the raw and the corrected, after the subtraction of the gamma and multiple scattering background, TOF NCS data for d-isopropanol and isopropanol, respectively. It is worth noting that upon subtraction of both types of background the magnitude of the errorbars of the TOF NCS data does not increase. The reason for this is twofold. Firstly, the gamma-background correction is a result of the application of an analytical, and thus smooth, function of the time of flight [46]. Secondly, the multiple scattering correction, being a result of a Monte-Carlo-type calculation, can be made arbitrarily noiseless by increasing the number of events [24]. In practice, already one million individual neutron histories propagated through the Monte-Carlo simulation produces practically noiseless signals [24] . Moreover, the errors of the fitting parameters (peak widths and intensities) almost do not propagate through the gamma-background and multiple-scattering corrections, largely owing to the fact that they result from integration-type of calculations [24,46]. In practice, the shapes of both types of corrections do not change within the two standard deviations of the uncertainties of the fitting parameters, which is tantamount to the iterative approach to TOF NCS data correction and fitting [24,46].
The coloured solid lines in the top and middle panel of figure 2 represent the fitting curves of the fitting curves calculated using the expression for the total count rate for TOF NCS data (given by Eg. (1)), as well as the fitting curves representing contributions of recoil peaks of individual masses to the total count rate. As mentioned in point 2 of the TOF NCS data treatment protocol described above, the fitting curves were computed using the stoichiometrically fixed relative peak intensities. The technique of stoichiometric fixing is a potent tool in fitting TOF NCS data as it enhances the precision of the determination of the recoil peak widths being the only non-linear fitting parameters of the model given by Eg. (1). Importantly in the context of the determination of the NMD widths and kinetic energy values of deuteron, the stoichiometric fixing greatly improves the precision of the determination of the width of the deuteron NMD. When averaged over results obtained from all forward scattering detectors, the uncertainty of the deuteron NMD width does not exceed 0.2 Å −1 , and this value will be used hereinafter as an upper conservative bound. The bottom panel in figure 2 shows the longitudinal momentum distribution of the hydrogen in isopropanol, the hydrogen in d-isopropanol, and their difference, representing the longitudinal momentum distribution of the hydrogen in the OH group, together with the fitting curve and the hydrogen resolution function in the longitudinal momentum space. The longitudinal momentum distribution of the hydrogen has been scaled by the factor of 7/8 in order to account for a different number of protons in isopropanol and d-isopropanol. It is worth noting that this scaling can only be performed in the longitudinal momentum space of the proton, not in the TOF domain. The reason for this is that the TOF NCS spectra contain only the information about the relative number of moles of the constituent nuclear species per formula unit of the sample (see equation (1)). Thus, a subtraction of a normalised and scaled TOF NCS spectrum of d-isopropanol from the normalised TOF NCS spectrum of isopropanol would inevitably lead to a creation of a negative residual signal in the TOF region of the deuteron recoil peak and a positive residual signal in the TOF region of the carbon, oxygen, and aluminium recoil peaks. In other words, one cannot apply a single scaling factor to the entire TOF NCS spectrum in order to isolate the TOF NCS contribution from the hydrogen in the OH group of the isopropanol. Such operation must be performed after the recoil peak isolation (subtraction of the recoil peaks of all other nuclear species apart from the hydrogen) in the TOF domain. Moreover, the isolated in the TOF domain proton recoil peaks in isopropanol and d-isopropanol correspond to different scattering angles and thus are not centred at the same TOF values, thus rendering the peak subtraction less precise an operation. Thus, focusing on the proton longitudinal momentum space followed by peak subtraction in the same space seems like the most numerically stable and precise option.
Finally, let us mention that the longitudinal momentum distribution of the OH group proton, illustrated as a solid blue line in the bottom panel in figure 2, can be fitted in a more reliable, robust, and precise manner owing to the signal differencing procedure outlined above. Namely, all signal fluctuations other than random Poissonian noise that is present in both longitudinal momentum distributions of isopropanol and d-isopropanol cancel out in this differencing procedure. This is clearly seen in the data shown in the region of longitudinal momentum magnitudes between 5 and 10 inverse Angstroms, on the right-hand side of the bottom panel in figure 2. The resultant difference-curve is then much smoother and more regular, resulting in more precise values of the fitted parameters. Moreover, the resolution function for the measurement of the proton longitudinal momentum distribution (red solid area in the bottom panel of figure 2) is much narrower than in the case of other nuclear mass-dependent VESUVIO resolution functions [8, 10-12, 47, 48], thus enabling highprecision work. Taken together, data focusing and high-resolution longitudinal nuclear momentum distribution spectroscopy yield withs of proton momentum distributions in isopropanol and d-isopropanol with the precision better than 0.2 inverse Angstrom. Thus, a value of 0.2 inverse Angstrom will be hereafter used as an upper conservative bound.

Results and discussion
As a first step in the interpretation of the values of the recoil peak widths, obtained from the NCS experiments, an ab initio calculation was performed. Two sets of simple ab initio calculations with the Density Functional Theory (DFT) method B3LYP using the basis set B3LYP/6-31 g * using Gaussian 98, revision A11, were performed in the limit of a single isolated molecule: (i) in vacuo calculation, and (ii) a calculation in a solvent environment of water. In the latter case, the solvent was modelled by placing the molecule within the solvent reaction field. For this, the self-consistent reaction field (SCRF) methodology, implemented using the COnductor-like Screening MOdel (COSMO) variant [49] of the CPCM model [50], was applied. COSMO is a variant of the continuum solvation models, which uses a scaled conductor boundary condition instead of the much more complicated dielectric boundary condition for the calculation of the polarization charges of a molecule in a continuum. The optimisation procedure led to a stable configuration and no imaginary frequencies. Table 1 shows the predictions for the isotropic average values of the widths of the proton, carbon and oxygen momentum distributions in isopropanol based on the DFT calculation in vacuo and the calculation in a solvent environment of water. Table 2 shows the results of the analogical calculation for d-isopropanol, also including the value for D in the OD group. First, as an aside, it is worth noting that, in the case of the carbon, oxygen, and deuteron, the above theoretical predictions encourage the conclusion that, in principle, the NCS method is capable of achieving site and solvent selectivity in the idealised case of nearly-perfect resolution.
Taking into account, however, that the single-standard deviation (1-STD) error-levels of the present measurements of the widths of C and O momentum distributions are around 0.7-0.8 Å −1 , this remains still in the realm of the future vision for the much-needed upgrade of the NCS technique. Thus, site-selective isotopic substitution remains at present the only plausible route to achieving site-selectivity of the NCS method beyond the proton, a realisation that motivated this work.
The main observation from the DFT calculation mentioned above is that, unlike for carbon, for hydrogen, the DFT-predicted values, both in vacuo and in the presence of the water solvent model, are almost independent on the position of the hydrogen atom in the molecule. Treating the sets of these values as distributions, one can conclude that the standard deviations of these distributions amount to less than 1% of their mean values. This is an encouraging result, bearing in mind that NCS is not site-selectively measuring the proton momentum distributions. It means that in this particular case one can compare the widths of proton momentum distributions from NCS experiments with doubly-averaged (first over directions in space for each individual type of proton, and then over different types of protons in the molecule) theoretical predictions. Thus, one can assume that any discrepancy between the average theoretical prediction and the experimental value must signal a failure of the harmonic approximation to account for the local-effective potentials experienced by all protons in isopropanol and d-isopropanol. In what follows, our discussion about the degree of anharmonicity in both molecular systems will be motivated by this notion. Figure 3 shows the values of the widths of nuclear momentum distributions of H in isopropanol, H in d-isopropanol, and D in d-isopropanol, obtained from fitting Gaussian profile functions to recorded NCS data, as well as the theoretical predictions for the average H and D widths, based on the harmonic ab initio DFT in vacuo calculation and the calculation in a solvent environment of water.     figure 3 are systematically below their respective DFT predictions, calculated within the harmonic approximation (HA). Moreover, the viral theorem predicts, in the case of pronounced anharmonicity of the local effective potential experienced by the nuclei, that the average nuclear kinetic energy becomes less than half of the total energy of the particle confined in the potential [51]. Taken together, the results shown in figure 3 signal a failure of the harmonic approximation to account for the local-effective potentials experienced by all protons in isopropanol and d-isopropanol. The next step in our appraisal of the local effective potential experienced by H and D in isopropanol and d-isopropanol is summarised in figure 4. The top panel of figure 4 shows the ratios of the kinetic energies of H in OH in isopropanol and D in OD in d-isopropanol, obtained from the widths of NMDs using the expression: It is worth noting that, by virtue of error propagation, the error of the ratio of the kinetic energies given by Eg. It is worth noting that the theoretical ratio plotted in figure 4, Ekin(H in OH)/Ekin(D in OD), was not taken to be the ratio of the reduced masses of D in OD and H in OH, being 1.373. The reason for that is that the NCS method is probing the total nuclear kinetic energies not the energies of the individual modes of vibration of the nuclei. The total nuclear kinetic energies are functions of the Boltzmann-population-weighted first-moments of the atom-projected VDOS distributions. These quantities are in turn calculated taking the masses of the freely recoiling nuclei within the Impulse Approximation. Physically, the values of these masses are given by the conservation of the total kinetic energy and momentum of the system consisting of an impinging neutron and a target nucleus that recoils after the collision with the neutron during the Compton scattering event. Mathematically, the conservation of the nuclear kinetic energy and momentum during a Compton event is reflected in the fact that the longitudinal momentum distributions of all nuclear species are centred at the values of zero. This fact is clearly seen in the bottom panel of figure 2 for H in OH. The centre of the longitudinal momentum distribution for the proton in the OH group in this figure corresponds to the zero value of the longitudinal momentum as calculated for the value of the freely recoiling proton of the mass of 1.0079 amu. For the case of D in OD, the same property is reflected in the fact that the TOF NCS data treatment of the recoil peaks of D in d-isopropanol always results in the recoil peak positions in the TOF domain that correspond to the value of zero longitudinal nuclear momentum of the freely recoil mass of D, i. e., 2.015 amu (see the top panel in figure 3). The physical picture pertaining to the momentum distributions of H in OH and D in OD is thus that they both correspond to the distributions of all nuclear momentum components of all modes of vibrations that have non-zero contributions to the motions of H in OH and D in OD functional groups. Although the bulk of contributions to these distributions comes from the nuclear stretching modes, there are also non-negligible contributions from other modes of vibrations, such as bending modes for example. Thus, there is no single effective mode of vibration and a single values of effective mass pertaining to this mode that is representative for the nuclear motion of H in OH and D in OD. In other words, the effective mean forces acting on of H in OH and D in OD are not aligned along the directions of the OH and OD bond respectively but rather reflect complicated interplays between all vibrational modes. It is worth noting that, the theoretical predictions taking into account Boltzmann population factors yield kinetic energy isotope effect that is decreasing with increasing temperature. This is the reflection of the fact that, as temperature increases, the high-energy vibrational modes contributing to the motion of D in OD are being populated faster than the modes in the case of H in OH. In consequence, the temperature-weighted mean values of the atom-projected vibrational densities of states (and thus the average kinetic energy values) shift towards higher energies faster in the case of D in OD. In other words, the trends represented by solid black and blue lines, when extrapolated to T=0 K, would reach the value of 1.41. As an aside, we remark that an interesting solvent effect is clearly visible in theoretical prediction depicted in figure 3. Namely, the presence of the solvent (water) seems to shift the theoretical prediction towards lower magnitudes, as compared to the prediction in vacuo. This indicates that, unlike the absolute values of the nuclear kinetic energy for H and D, the site-selective values of the magnitude of the isotope effect for the kinetic energies can potentially be used as observables of choice for the detection of solvent-induced nuclear quantum effects. Of course, this observation is of preliminary character and needs further confirmations by means of more sophisticated ab initio methodology (e. g., ab initio molecular dynamics or simulations). Notwithstanding the last remark, a clear departure is visible in figure 4 for the NCS data from the harmonic approximation results. Firstly, the experimental values are placed systematically higher than any harmonic vibration prediction. Secondly, the temperature trend for the experimental data is opposite as compared to the theoretical predictions. This result further corroborates the notion of the anharmonic character of the local-effective potentials of H in OH in isopropanol and D in OD in d-isopropanol in the condensed phase.
Further insight into the possible shape of the local effective potential experienced by H in OH and D in OD in condensed phases of isopropanol and d-isopropanol can be gained owing to the analysis of trends visible is spectroscopic measurements for systems with hydrogen bonds, performed by McKenzie et al [52]. The first observation is that the magnitudes of the isotope effect for the nuclear kinetic energy, observed in our NCS experiments, are well within the range of values observed for the longitudinal O-H(D) stretching mode for hydrogen to deuteron isotopes is observed experimentally, varying from 0.85 to 2.0 [52][53][54][55]. In contrast, for the torsional/bending modes, the isotope effects are consistent with the semi-classical harmonic ratio of 1.41 [52][53][54][55]. Thus, the vibrational frequency isotope effect is expected to correlate well with the nuclear kinetic energy isotope effect for H(D) in hydrogen bonds like those formed by H(D) in OH (OD) in isopropanol and d-isopropanol. Both systems investigated here belong to the class of weakly hydrogen-bonded systems. Unlike the nearly harmonic weak hydrogen bonds, moderate to weak hydrogen bonds are anharmonic at the bottom of their effective BO potential energy curve. This type of anharmonicity can be very well parametrised using the Morse potential [52]. The Morse potential energy function is of the form: where r is the distance between the atoms, r e is the equilibrium bond distance, D is the potential well depth (defined relative to the dissociated atoms), and a is the attenuation length that controls the width of the potential. The Morse potential, rather than being defined by equation (6), is often rewritten in a form in which the zero of the potential energy is shifted down by subtracting the value of the dissociation constant D, which gives: The Eigen energies of the Schrodinger equation with the Morse potential can be written in terms of its original variables as: This procedure is shown in figure 5 for two temperatures, at which NCS were performed, T=250 K and T=300 K. Assuming a and D fulfil equation (9) at both T=250 K and T=300 K, we obtain a=2.2±0.2 Å −1 and D=6713±515 meV. Interestingly, the obtained values of a and D fall within the range of Morse potential parameters tabulated in the literature for simple molecular systems [56].
In the last step of our appraisal of the local-effective potential experienced by H in OH in isopropanol, we calculate the degree of the delocalisation in space and momentum of the muoniated variant of the isopropanol. We start by noting that it is plausible to assume that the Born-Oppenheimer approximation (BOA) is valid for the case of H, D, and Mu-substituted molecular species. The reasons to assume this are twofold. Firstly, the assumption of the validity of the BOA and the lack of non-adiabatic effects have already been tested both theoretically [57] and experimentally [58,59] for the case of the neutron Compton scattering on lightweight nuclei. Mathematically, the theory predicts that the BOA violation in NCS would be reflected in a build-up of satellite peaks at the low-energy transfer shoulders of the main recoil peaks (that are present in the ideal case of the absence of non-BOA effects) [57]. Such non-BOA effect would be most pronounced in the case of lightweight nuclei [57]. In the longitudinal momentum distribution domains, such non-BOA effects would have led to shifts of the centroids of nuclear longitudinal momentum distributions of H and D towards lower values of the longitudinal momentum [57]. However, no such effects have been experimentally observed in the NCS spectroscopy [58,59]. Such is the case also for the nuclear momentum distributions of H in isopropanol and D in d-isopropanol discussed in this work. Namely, no additional shifts of the centroids towards lower values of the longitudinal nuclear momentum have been observed here during the NCS data treatment described above.
When the hydrogen is substituted by a positive muon, the nucleus-to-electron mass ratio falls from a value of 1836.15 to 206.77 with the resulting atom muonium, Mu, behaving like a light isotope of hydrogen [60]. One could anticipate that for muonium-substituted molecules, the isotope effect is so large that the BOA is invalid. However, non-adiabatic and non-BOA effects are extremely small as observed by the muon-spin resonance and predicted by theoretical calculations for Mu-substituted molecular species and solid-state defect centres, at least in their ground electronic states [61]. For example, for isotopomers of the hydrogen molecule Webster and McKenna showed, by a series of non-adiabatic calculations beyond the Born-Oppenheimer approximation, that the key isotope effect arises through the mass dependence of the zero-point vibrational energy but essentially the same single potential-energy curve underlies the nuclear dynamics of all nuclear isotopic species [62]. On the whole, this result indicates the so-called trivial isotope effect, which, in order to be accounted for, in essence does not require resorting to non-adiabatic calculations. The reason for this is that the electronic structure for the protonic and muonic species is identical within the Born-Oppenheimer approximation, mass effects being manifest only through differences in zero-point vibrational motion and the alteration of the vibrational modes. A further theoretical study of the muonium-substituted species MuCO · and MuCO + also points out on the trivial isotope effect being the mere consequence of the mass ratio of the isotopic species. It was namely indicated that the vibrational frequencies of the C-Mu bond in the radical MuCO · increase by a factor close to a value of three, as expected for harmonic vibrations where the frequency varies inversely with the square root of the reduced mass for a Mu-C bond [63]. Furthermore, another study of vibrational behaviour of the carbon-muonium bond in the muoniated ethyl radical points out that, despite the great disparity in mass of the hydrogen and muonium isotopomers, the mere effect of the isotopic substitution is that the zero-point vibrational level lies significantly higher on the potential surface for Mu and the C-Mu bond executes large amplitude vibrations compared with a C-H bond, with the mean-square displacement for the C-Mu bond being three-time larger than the figure for the C-H bond [60]. However, using the relation between the magnitude of the force constant acting on the nuclear species in a bond, the magnitude of the mean-square displacement and the value of the zero-point energy for the stretching mode of the C-H bond, the authors arrive at the conclusion that values of the equilibrium stretching force constant are essentially isotope independent, despite the fact that the C-Mu bond executes large-amplitude vibrations and the zero-point vibrational energy of C-Mu is large compared with a C-H bond [60]. Using the language of the NCS one could thus state that the effective mean forces acting on Mu in the C-Mu bond and H in the C-H bond and their respective potentials are essentially the same. This is a significant observation as it essentially lies at the foundation of the so-called steric isotope effect [60]. Namely, the isotopic substitution introduces overcrowding if, as is supposed, a muonium atom has a larger steric requirement than hydrogen [60]. The energetic consequence of the steric isotope effect is that there is an isotopic shift, induced by steric interaction, in the zero-point vibrational energy [60].
Taken together, the experimental and theoretical results discussed above point out that one can assume that muon and proton are governed by almost identical potential-energy surfaces in analogous species [61]. Moreover, the notion of the equilibrium or Born-Oppenheimer structure can be used in the above considerations as a benchmark for the estimate of the magnitude of the isotope effect as long as the BOA is considered valid in the muoniated radicals. Importantly, the BOA leads to a potential energy surface that is a function of nuclear positions but not of nuclear masses. The validity of the BOA for Mu-substituted molecules has often been questioned on the qualitative level. However, as we have seen in the discussion in the last few paragraphs, the quantitative theoretical assessments lead to the conclusion that the departure from the BOA manifests as an increase in force constant confining the muon of not more than a few per cent. The effect on observable molecular properties, such as the vibrational frequencies, root-mean-square displacements and the widths of nuclear momentum distributions, varying as / k m , -( ) / mk 1 4 and ( ) / mk 1 4 respectively, are likely to be very small [61]. Thus, the isotope effects are primarily due to zero-point energy in the ground vibrational state.
Assuming that both hydrogen and Mu in OH in isopropanol experience the same Morse potential, one can calculate the kinetic energies of both particles according to equation (9). Taking the parameters of the potential as a=2. =0.34, thus a value higher than one third, assumed in the case of the harmonic approximation and also a result typical for the Morse potential [61]. Using the definition of the kinetic energy (equation (3)) and assuming that the width of the proton momentum distribution in OH is 4.7 Å −1 , we get the prediction for the width of the Mu momentum distribution of 2.7 Å −1 . Inverting between the width of the momentum and spatial delocalisation (using the Heisenberg uncertainty principle for Gaussian wave packets), we get the prediction for the extent of the spatial delocalisation of the Mu in the muoniated isopropyl radical of 0.37 Å.
It is tempting to put our estimates of the kinetic energy, space and momentum delocalisation of the Mu in isopropyl alcohol in the context of the literature results on nuclear quantum and spin dynamics of Mu in condensed matter. Mu, when embedded in solid-state systems is usually shown to be dominated by the tendency to form bonds with adjacent atoms. For example, for Mu in silicon the values of the zero-point energy are, depending on its position between the two Si atoms in the lattice, between 240 and 270 meV, for the motion along the bond direction, and between 540 and 560 meV when the motion perpendicular to the bond is included [64]. These amounts of zero-point energy are sufficient to place the average position of the muon at the siliconsilicon bond centre. However, because of the extraordinary light mass of the muon, the inclusion of the zeropoint motion expands the bond by 0.88 Å. Moreover, if one replaces the muon with a proton the corresponding extension is about 0.86, and 0.84 Å for the deuteron [64]. From our estimates we are unable to predict the magnitude of the hydrogen bond extension in the muoniated isopropyl. We can, however, assess the amount of the zero-point energy of the Mu. Taking the kinetic energy of the proton in OH in isopropanol as 132.5 meV at T=300 K and scaling it by 1/0.34 we obtain an estimate of 389.8 meV for the kinetic energy and 779.6 meV for the zero-point energy of the Mu in the muoniated isopropanol at the same temperature. This value is ca. 40% larger than the value reported for the Mu in the silicon lattice. Thus, one can anticipate that this amount of the zero-point energy may elongate the hydrogen bond in the muoniated isopropyl beyond the value of 0.88 Å.
The magnitude of the zero-point energy can also be used to predict the transition from the classical to the quantum diffusion regime of Mu in solid-state systems. It is predicted that the transition between the classical and quantum Mu diffusion takes place at a temperature T * ∼ν 0 , where ν 0 is the vibrational frequency of the particle in its potential well [65]. For example, for Mu atoms in crystals, T * can be estimated to be on the order of several hundred Kelvin [66,67]. Inspired by this result, we can provide an upper conservative bound for the temperature corresponding to the onset of the classical hopping of the Mu in the isopropyl in the liquid phase at ambient conditions. Our estimate that the magnitude of the zero-point energy of the Mu in the liquid isopropanol is of the order of ca. 9000 K and 188.4 THz, respectively, 2 orders of magnitude more than the values reported for solids. Thus, quantum diffusion of the Mu is much more likely at ambient conditions in the liquid isopropanol than in many solid-state systems.
Last but not least, the isotope effects on the zero-point energy between Mu, hydrogen, and deuteron in molecular and solid-state systems are widely discussed in the context of the validity of the predictions for the magnitudes of the hyperfine coupling constants in μSR work [61], a discussion that originally motivated the NCS work presented here. As a starting remark it is worth noticing that the classical or equilibrium structures of different isotopomers of the same system would be the same in the Born-Oppenheimer approximation, leading to identical equilibrium bond-lengths and identical spin density on deuteron, proton or muon, with their hyperfine constants scaling according to their gyromagnetic ratios [61]. However, experimentally this is not the case, both for rigid or more flexible molecules, a natural consequence of the different zero-point energies for muon and proton, and deuteron [61]. The relatively rigid structures of C60X (the monohydride isotopomers with X=D, H, Mu) provide a simple introduction to hyperfine isotope effects [61]. The influence of the different zero-point energies of the three isotopes on the magnitude of the hyperfine coupling is even more pronounced if the local-effective potential energy curve is anharmonic, such as the one describing D, H, Mu in C60X, which is very well approximated by the Morse potential [61]. Contrary to the harmonic case, there are now two main contributions to the isotope effect on the hyperfine coupling, the harmonic and the anharmonic one. Stretching and rocking modes as well as the motion transverse to the original double bond for C60X project on the motion of H and Mu differently. The net effect is that the position distribution of Mu is much more asymmetric and skewed towards the softer side of the potential than the position distribution of H. In such a scenario, the anharmonic contribution comes from the different average bond lengths and the harmonic contribution from the different root-mean-square excursions of H and Mu from the minimum of the potential energy surface [61]. In consequence, the different amounts of zero-point energy for H and Mu, equal to 220 and 590 meV respectively, are distributed in space completely differently, leading to completely different sampling by H and Mu of the unpaired electron spin density. One can imagine that, as the Mu bond is stretched to breaking point, hyperfine constants increase to the respective atomic values [61]. On the whole, the spin density on the deuteron, proton and muon is higher than the equilibrium-geometry value by 2%, 3% and 10%, respectively, leading to overall isotope effects of 1% between C60D and C60H, and about 7% between C60H and C60Mu. Taking into account that our estimate of the zero-point energy of Mu in isopropyl is ca. 32% larger than the value for C60Mu, we can anticipate that the overall isotope effect for the hyperfine interaction between isopropyl and its muoniated counterpart is expected to be at least one third larger, i.e. in the region of 10%. Moreover, one can treat this estimate as lower conservative bound, as it assumes that the isopropyl is a rigid molecular system. However, in muoniated organic radicals with softer internal modes such as hindered rotations, hyperfine isotope effects can reach up to 20%. These large magnitudes of the isotope effects are mainly due to the averaging over bond angles (largely absent in the case of rigid molecules) and the role played by the soft modes [61].
We close our discussion by remarking that the estimates of the kinetic and zero-point energy, as well as the position and momentum delocalisation of the Mu in the muoniated isopropyl molecule in condensed phase is by no means an exact calculation, but merely a first step towards further elucidation of the nuclear and spin dynamics of Mu in this important system. One can imagine that, in the future, much more precise ab initio work, based for example on the PIMD approach, will be available. To this end, one can envisage a computational protocol in which our estimates are treated as additional constraints for the PIMD calculation, hopefully guiding it towards a much better agreement with experiment.

Summary
In this work, we have employed the combination of NCS and ab initio single-molecule vibrational calculation in order to perform a critical appraisal of the nuclear quantum dynamics of protons and deuterons in isopropyl and d-isopropyl in the condensed phase. Based on the classical results for the magnitude of the isotope effect for model inter-particle potentials, we were able to assess the degree of the isotope effect on nuclear and spin dynamics of Mu in muoniated isopropyl. The NCS technique has been employed for the first time as a tool shedding more light on the nuclear and spin dynamics of the Mu in the condensed phase. The obtained results can serve as conservative constraints for future precise ab initio calculations of properties of muoniated radicals in condensed matter systems and molecules.