Transfer matrix theory of surface spin echo experiments with molecules

TL;DR

This paper develops a quantum mechanical transfer matrix method to model molecular surface spin-echo experiments, enabling the extraction of surface scattering information from experimental signals involving molecules like ortho-hydrogen.

Contribution

It introduces a one-dimensional transfer matrix approach that incorporates molecular hyperfine states and magnetic field effects, advancing the theoretical modeling of molecular surface spin-echo experiments.

Findings

The method accurately connects experimental signals to scattering matrix elements.

Preliminary comparison shows sensitivity of signals to surface scattering properties.

Framework enables future Bayesian inference of surface dynamics.

Abstract

${}^3$He beam spin-echo experiments have been used to study surface morphology, molecular and atomic surface diffusion, phonon dispersions, phason dispersions and phase transitions of ionic liquids. However, the interactions between ${}^3$He atoms and surfaces or their adsorbates are typically isotropic and weak. To overcome these limitations, one can use molecules instead of ${}^3$He in surface spin-echo experiments. The molecular degrees of freedom, such as rotation, may be exploited to provide additional insight into surfaces and the behaviour of their adsorbates. Indeed, a recent experiment has shown that $\textit{ortho}$-hydrogen can be used as a probe that is sensitive to the orientation of a Cu(115) surface [Godsi et al., Nat. Comm. $\textbf{8}$, 15357 (2017)]. However, the additional degrees of freedom offered by molecules also pose a theoretical challenge: a large manifold of molecular states and magnetic field-induced couplings between internal states. Here, we present a fully quantum mechanical approach to model molecular surface spin-echo experiments and connect the experimental signal to the elements of the time-independent molecule-surface scattering matrix. We present a one-dimensional transfer matrix method that includes the molecular hyperfine degrees of freedom and accounts for the spatial separation of the molecular wavepackets due to the magnetic control fields. We apply the method to the case of $\textit{ortho}$-hydrogen, show that the calculated experimental signal is sensitive to the scattering matrix elements, and perform a preliminary comparison to experiment. This work sets the stage for Bayesian optimization to determine the scattering matrix elements from experimental measurements and for a framework that describes molecular surface spin-echo experiments to study dynamic surfaces.

Tables

Table 1. Table 1: Relative probabilities of the state selector η m I ​ m J subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 \eta_{m_{I}m_{J}} and the detector κ m I ​ m J subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 \kappa_{m_{I}m_{J}} . The state selector probabilities P R 0 subscript 𝑃 subscript 𝑅 0 P_{R_{0}} are calculated as P R 0 = P m I ​ m J ≡ η m I ​ m J / ∑ m I ​ m J η m I ​ m J subscript 𝑃 subscript 𝑅 0 subscript 𝑃 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 P_{R_{0}}=P_{m_{I}m_{J}}\equiv\eta_{m_{I}m_{J}}/\sum_{m_{I}m_{J}}\eta_{m_{I}m_{J}} and the detector coefficients c R D subscript 𝑐 subscript 𝑅 𝐷 c_{R_{D}} are calculated as c R D = c m I ​ m J ≡ κ m I ​ m J / ∑ m I ​ m J κ m I ​ m J subscript 𝑐 subscript 𝑅 𝐷 subscript 𝑐 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 c_{R_{D}}=c_{m_{I}m_{J}}\equiv\kappa_{m_{I}m_{J}}/\sum_{m_{I}m_{J}}\kappa_{m_{I}m_{J}} .

$m_I$	1	1	1	0	0	0	-1	-1	-1
$m_J$	1	0	-1	1	0	-1	1	0	-1
${\eta }_{m_I{m}_J}$	0.0095	0.0138	0.0187	0.0416	0.0436	0.0606	0.3997	0.9015	1.0
${\kappa }_{m_I{m}_J}$	0.0611	0.08	0.1027	0.3834	0.5705	0.8425	1.0	0.9422	0.7209

Table 2. Table 2: Relative probabilities of the state selector η m I ​ m J subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 \eta_{m_{I}m_{J}} and the detector κ m I ​ m J subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 \kappa_{m_{I}m_{J}} , as used in the comparison to the semi-classical method of Godsi et al.  Godsi  et al. ( 2017 ) . The state selector probabilities P R 0 subscript 𝑃 subscript 𝑅 0 P_{R_{0}} are calculated as P R 0 = P m I ​ m J ≡ η m I ​ m J / ∑ m I ​ m J η m I ​ m J subscript 𝑃 subscript 𝑅 0 subscript 𝑃 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜂 subscript 𝑚 𝐼 subscript 𝑚 𝐽 P_{R_{0}}=P_{m_{I}m_{J}}\equiv\eta_{m_{I}m_{J}}/\sum_{m_{I}m_{J}}\eta_{m_{I}m_{J}} and the detector coefficients c R D subscript 𝑐 subscript 𝑅 𝐷 c_{R_{D}} are calculated as c R D = c m I ​ m J ≡ κ m I ​ m J / ∑ m I ​ m J κ m I ​ m J subscript 𝑐 subscript 𝑅 𝐷 subscript 𝑐 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript subscript 𝑚 𝐼 subscript 𝑚 𝐽 subscript 𝜅 subscript 𝑚 𝐼 subscript 𝑚 𝐽 c_{R_{D}}=c_{m_{I}m_{J}}\equiv\kappa_{m_{I}m_{J}}/\sum_{m_{I}m_{J}}\kappa_{m_{I}m_{J}} .

$m_I$	1	1	1	0	0	0	-1	-1	-1
$m_J$	1	0	-1	1	0	-1	1	0	-1
${\eta }_{m_I{m}_J}$	1.0000	0.9755	0.7901	0.1465	0.1111	0.0738	0.0343	0.0299	0.0258
${\kappa }_{m_I{m}_J}$	1.00	0.96	0.93	0.53	0.42	0.37	0.21	0.19	0.16

Equations

$${\hat{H}}^{\text{R}}(\vec{B})=U\left({\hat{I}}^2,{\hat{J}}^2,\hat{I}\cdot \hat{J},I,J\right)+V\left({\hat{I}}^2,\hat{I}\cdot \vec{B},I,{\vec{B}}^2\right)+Q\left({\hat{J}}^2,\hat{J}\cdot \vec{B},J,{\vec{B}}^2\right),$$

$$\hat{H}(x)=\frac{{\hat{p}}^2}{2m}+{\hat{H}}^{\text{R}}(\vec{B}(x))$$

$${\rho }_0$$

$$|{\Psi }_{R_0{k}_0}⟩$$

$${\psi }_{R_0{k}_0}(r)$$

$$=\sqrt{{\sigma }_k}{\left(\frac{2}{\pi }\right)}^{\frac{1}{4}}{e}^{i{k}_0^{R_0}r}{e}^{-r^2{\sigma }_k^2}$$

$$P_{\text{detection}}=\frac{1}{\tau }{\int }_{t_1}^{t_2}\text{d}t⟨\hat{C}(t)⟩,$$

$$⟨\hat{C}(t)⟩=\mathrm{Tr}\hat{\rho }(t)\hat{C},$$

$$\hat{C}=\sum _{R_D}\int \text{d}x\text{leavevmode}{c}_{R_D}(x)|x{R}_D⟩⟨x{R}_D|$$

$$\rho (t)=\sum _{R_{0}R{R}^{\prime }}\int \text{d}E\text{leavevmode}\int \text{d}{E}^{\prime }\text{leavevmode}{P}_{R_0}{e}^{-\frac{i}{\hslash }(E-E^{\prime })t}{\alpha }_{k_0{R}_0}^{ER}{\alpha }_{k_0{R}_0}^{\ast E^{\prime }{R}^{\prime }}|ER⟩⟨E^{\prime }{R}^{\prime }|,$$

$$⟨\hat{C}(t)⟩=\sum _{R_D,R_D^{\prime }}\int \text{d}r\text{leavevmode}\int \text{d}{r}^{\prime }\text{leavevmode}⟨r^{\prime }{R}_D^{\prime }|\rho (t)|r{R}_D⟩⟨r{R}_D|\hat{C}|r^{\prime }{R}_D^{\prime }⟩$$

$$⟨r^{\prime }{R}_D^{\prime }|\rho (t)|r{R}_D⟩=\sum _{R_{0}R{R}^{\prime }}\int \text{d}E\text{leavevmode}\int \text{d}{E}^{\prime }\text{leavevmode}{P}_{R_0}{e}^{-\frac{i}{\hslash }(E-E^{\prime })t}{\alpha }_{k_0{R}_0}^{ER}{\alpha }_{k_0{R}_0}^{\ast E^{\prime }{R}^{\prime }}{\Phi }_{R_D^{\prime }}^{ER}(r^{\prime }){\Phi }_{R_D}^{\ast E^{\prime }{R}^{\prime }}(r).$$

$$⟨r{R}_D|\hat{C}|r^{\prime }{R}_D^{\prime }⟩$$

$$=c_{R_D}(r)\delta (r^{\prime }-r){\delta }_{R_D^{\prime }{R}_D}$$

$$⟨\hat{C}(t)⟩=\sum _{R_{0}R{R}^{\prime }}\int \text{d}E\text{leavevmode}\int \text{d}{E}^{\prime }\text{leavevmode}{P}_{R_0}{e}^{-\frac{i}{\hslash }(E-E^{\prime })t}{\alpha }_{k_0{R}_0}^{ER}{\alpha }_{k_0{R}_0}^{\ast E^{\prime }{R}^{\prime }}\left(\sum _{R_D}\int \text{d}r\text{leavevmode}{\Phi }_{R_D}^{ER}(r){\Phi }_{R_D}^{\ast E^{\prime }{R}^{\prime }}(r)c_{R_D}(r)\right).$$

$${\alpha }_{k_0{R}_0}^{ER}$$

$$={\delta }_{R{R}_0}{\Gamma }_{k_0{R}_0}^{ER}$$

$$⟨\hat{C}(t)⟩$$

$$P_{\text{detection}}$$

$$=\sum _{R_0}\int \text{d}E\text{leavevmode}\int \text{d}{E}^{\prime }\text{leavevmode}{P}_{R_0}\frac{2\pi \hslash }{\tau }\delta \left(E-E^{\prime }\right){\Gamma }_{k_0{R}_0}^{E{R}_0}{\Gamma }_{k_0{R}_0}^{\ast E^{\prime }{R}_0}\left(\sum _{R_D}\int \text{d}r\text{leavevmode}{\Phi }_{R_D}^{E{R}_0}(r){\Phi }_{R_D}^{\ast E^{\prime }{R}_0}(r)c_{R_D}(r)\right)$$

$$=\sum _{R_0}\int \text{d}E\text{leavevmode}{P}_{R_0}{\left|{\Gamma }_{k_0{R}_0}^{E{R}_0}\right|}^2\left(\sum _{R_D}\int \text{d}r\text{leavevmode}\frac{2\pi \hslash }{\tau }{\left|{\Phi }_{R_D}^{E{R}_0}(r)\right|}^2{c}_{R_D}(r)\right),$$

$${\Phi }_{R_D}^{E{R}_0}(r)$$

$$=e^{i{k}_{R_D}r}⟨R_D⟩E{R}_0$$

$$\equiv e^{i{k}_{R_D}r}{\beta }_{R_D}^{E{R}_0}$$

$$P_{\text{detection}}$$

$$|E\tilde{R}⟩=\sum _R\int \mathrm{d}x\text{leavevmode}{\Phi }_R^{E\tilde{R}}(x)|xR⟩,$$

$$-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}{\Phi }_{R_0}^{E\tilde{R}}(x)$$

$$\frac{{\partial }^2}{\partial x^2}{\Phi }_R^{E\tilde{R}}(x)$$

$${\Phi }_R^{E\tilde{R}}(x)=A_R{e}^{i{k}_{R}x}+B_R{e}^{-i{k}_{R}x},$$

$$k_R\equiv \frac{\sqrt{2m\left(E-E_{\text{R}}\right)}}{\hslash }.$$

$${\vec{\varphi }}_{x_2}={\mathbf{\Pi }}_{x_2-x_1}{\vec{\varphi }}_{x_1},$$

$${\mathbf{\Pi }}_x$$

$$=\left(\begin{array}{cccccc}{e}^{i{k}_{1}x}& & & & & \\ & \ddots & & & \text{makebox}(0.0,0.0)[]\text{0}& \\ & & e^{i{k}_{N_R}x}& & & \\ & & & e^{-i{k}_{1}x}& & \\ & \text{makebox}(0.0,0.0)[]\text{0}& & \ddots & & \\ & & & & & e^{-i{k}_{N_R}x}\end{array}\right),$$

Full text

I Introduction

A major thrust of recent experimental work has been to achieve control over the longitudinal motion of atomic and molecular beams Krems (2018); Lemeshko et al. (2013); Friedrich and Doyle (2009); van de Meerakker et al. (2008, 2012); Bethlem and Meijer (2003); Dulieu and Gabbanini (2009); Hogan et al. (2011). Controlled beams can be used for a variety of applications, ranging from loading molecules into traps Hutzler et al. (2012); Bethlem et al. (2000); Liu et al. (2015), to measuring cross sections for molecular scattering with extremely high energy resolution Onvlee et al. (2016); Jankunas and Osterwalder (2015); Stuhl et al. (2014); Brouard et al. (2014), to precision spectroscopy Tarbutt et al. (2013); Cheng et al. (2016); ACME Collaboration (2018), to controlled chemistry Segev et al. (2019). The development of methods for the initial state selection and control over both the longitudinal and transverse motion of molecular beams has also paved the way for matter-wave interferometry Cronin et al. (2009); Juffmann et al. (2013); Hornberger et al. (2012), nano-lithography Gordon et al. (2003); Rohwedder (2007); Dey et al. (2000) and precision studies of molecule-surface scattering. Although molecule-surface collisions have been a subject of numerous studies Benedek and Toennies (2018); Benedek and Valbusa (1982); Hulpke (1992); Tesa-Serrate et al. (2016), combining the latest advances in molecular beam control with surface scattering experiments opens opportunities for probing new regimes of molecule-surface energy exchange and obtaining detailed information about surface properties. This is well exemplified by 3He spin echo (HeSE) experiments DeKieviet et al. (1995, 1997); Jardine et al. (2009a, b) aiming to probe the structure of surfaces, as well as quantum matter adsorbed on surfaces, by scattering a beam of 3He in superpositions of nuclear spin states off a surface and observing the perturbation of the resulting interferometry signal. Analogous to neutron spin echo experiments Mezei (1980); Mezei et al. (2003), HeSE experiments have been shown to detect the impact of gravity (on the energy scale of $\approx$ 10 neV) on the kinetic energy of atoms in the beam DeKieviet et al. (1995). When used to study surfaces, HeSE experiments can be classified as a subset of quasi-elastic helium atom scattering experiments Benedek and Toennies (2018). Surface-sensitive HeSE experiments DeKieviet et al. (1997); Jardine et al. (2009a, b) have been used to study surface morphology Corem et al. (2013), molecular and atomic surface diffusion Jardine et al. (2009a, b); Hedgeland et al. (2016); Godsi et al. (2015); Hedgeland et al. (2011); Lechner et al. (2013a); Rotter et al. (2016), inter-adsorbate forces Godsi et al. (2015); Kole et al. (2012), phonon dispersions Kole et al. (2010); Jardine et al. (2009a, b), phason dispersions McIntosh et al. (2013), structures and phase transitions of ionic liquids McIntosh et al. (2014) and friction between adsorbates and surfaces Hedgeland et al. (2009a, b); Lechner et al. (2013b). HeSE experiments have provided information about potential energy surfaces Jardine et al. (2004, 2009a, 2009b) and surface-adsorbate interactions Jardine et al. (2009a, b, 2010) and are frequently combined with microscopic calculations to both test theory and gain insight into surface-adsorbate interactions Tamtögl et al. (2018); Sacchi et al. (2017); Hedgeland et al. (2011).

The use of 3He as probe particles in HeSE experiments can sometimes be limited by the weak interaction strength between 3He and surfaces or their adsorbates. In addition, 3He offers no internal degrees of freedom to absorb energy or induce anisotropic interactions. Therefore, an important recent goal has been to extend surface spin-echo experiments to molecular beams Godsi et al. (2017). Molecules offer rotational degrees of freedom and anisotropic, state-dependent interactions, which could be exploited to gain new insights into surface dynamics. For example, it was recently shown that ortho-hydrogen (oH2) can be used as a sensitive probe of surface morphology Godsi et al. (2017): the experiment was able to discern how the interaction between an oH2 molecule and a Cu(115) surface depends on the orientation of the rotational plane of the hydrogen molecule relative to the surface. In addition, one could exploit the transfer of rotational energy from the probe molecules to surface adsorbates (or vice versa) in order to study the relative effects of the rotational and translational motion on the dynamics of the adsorbates. However, the increased complexity of molecules (compared to 3He atoms) makes the analysis of the spin-echo experiments complicated and requires one to account for the interplay of the translational, nuclear spin and rotational degrees of freedom in strong magnetic fields of differing orientations, in addition to the molecule-surface scattering event.

Surface spin-echo experiments with molecules involve passing a molecular beam through a series of magnetic fields to control molecular wavepackets before and after the scattering event. A proper analysis of the resulting experimental signal must be based on (i) the solutions of the time-dependent Schrödinger equation accounting for the development of entanglement between the translational motion and the internal states of molecules in the beam, as the beam transverses the magnetic fields of the spin-echo apparatus; (ii) the description of the molecule-surface scattering events in the relevant frame of reference by the scattering matrix involving all relevant molecular states. This is a challenging task because the potential energy surfaces for molecule-surface interactions are difficult to compute with sufficient accuracy Golibrzuch et al. (2014); Krüger et al. (2015); Park et al. (2019); Yin et al. (2019); Jiang and Guo (2019); Tchakoua et al. (2019); del Cueto et al. (2019), the calculations of the cross sections for molecule-surface scattering are extremely time consuming Kroes et al. (2009); Kroes and Díaz (2016) and because the orientation and strength of magnetic fields necessarily change throughout the spin-echo apparatus. An alternative formulation can be developed to treat the molecule-surface scattering matrix elements as varying parameters to be determined from the experimental interferometry signal by one of the algorithms used in optimal control theory Gordon and Rice (1997); Balint-Kurti et al. (2008); Rabitz et al. (2000); Brif et al. (2010); Hofer and Hammerer (2017); Sola et al. (2018) or reinforcement machine learning designed to solve the inverse problem Vargas-Hernández et al. (2019); Rizzi et al. (2012). In order for such a formulation to be practical, it is necessary to develop a rigorous method for the description of molecular dynamics inside the spin-echo apparatus, before and after the molecular wavepackets interact with the surface. This method must be efficient to allow for multiple feedback control loops, be accurate to ensure the proper description of interferometry dynamics and integrate rigorously the surface scattering matrix amplitudes into the resulting output signal.

In this paper, we exploit the transfer matrix method Walker and Gathright (1992); Sánchez-Soto et al. (2012) to develop such a theoretical framework. The transfer matrix method Walker and Gathright (1992); Sánchez-Soto et al. (2012) has been applied in various fields, such as for solving the 2D Ising model in statisical mechanics Baxter (2007), calculating reflection and transmission coefficients in optics Born and Wolf (1980) and mesoscopic quantum transport Mello and Kumar (2004), determining photonic bandstructures Pendry and MacKinnon (1992), and examining the tunnelling of a molecule through potential barriers Saito and Kayanuma (1994); Jarvis and Bulte (1998). The general and efficient framework we present can be used to analyze the coherent propagation of closed shell molecules through a series of static magnetic fields with different magnitudes and orientations, as well as through one or more scattering events.

We apply this framework to surface-sensitive interferometry experiments that use closed shell molecules to study static surfaces. Specifically, we develop a fully quantum mechanical model of surface-sensitive molecular hyperfine interferometry experiments by deriving a one-dimensional transfer matrix method that includes the internal hyperfine degrees of freedom of the probe molecules and that accounts for the eigenbasis changes between local regions of the magnetic field. We account for the experimental geometry with rotation matrices and describe the molecule-surface interaction with a scattering transfer matrix (a transformed version of the standard scattering matrix).

The method is applied to an oH2 hyperfine interferometry experiment. By comparing the theoretical results with experimental measurements, we illustrate the importance of integrating over the velocity distribution of molecules in the beam. We further show that information about the scattering matrix elements is encoded in the experimental signal. In particular, we demonstrate that the experimental signal is sensitive to the magnitude and phase of the diagonal elements of the scattering transfer matrix. We also show that the signal is sensitive to scattering events that change the projection quantum numbers of the molecular hyperfine states. Such dynamical processes are described by scattering transfer matrices with non-zero diagonal and off-diagonal matrix elements. This sets the stage for determining, in part or in whole, the scattering transfer matrix elements of a particular molecule-surface interaction by comparing the computed and experimentally-measured signals.

Finally, we compare our method with a semi-classical method, which is described briefly in the supplementary material of Ref. Godsi et al. (2017) for oH2 and in more detail in Ref. Litvin et al. (2019) for spin 1/2 particles. Within this semi-classical method, the internal molecular degrees of freedom are treated quantum mechanically, while the centre of mass degree of freedom is treated classically. Through this comparison, we demonstrate that the present method can be extended to study dynamic, instead of static, surfaces by surface spin-echo experiments with molecules.

The remainder of this manuscript is organized as follows. In Section II, we describe a generic molecular hyperfine interferometry experiment. We then discuss, in Section III, the molecular state after the state-selecting magnetic lens. In Section IV, we time-evolve the molecular state and integrate the result over the length of the detection window to obtain the relationship between the system eigenstates and the detector current. To obtain the system eigenstates, we derive and apply, in Section V, a transfer matrix formalism that includes internal degrees of freedom. We also discuss the rotation and scattering transfer matrices used to account for the apparatus geometry and the molecule-surface interaction, respectively. In Section VI, we demonstrate the application of this theoretical framework to the case of oH2, illustrate the need to integrate over the velocity distribution, illustrate the sensitivity of the calculated signal to various features of the scattering transfer matrix and perform a preliminary comparison with experiment. We compare the method of the present manuscript to the semi-classical method discussed by Godsi et al. Godsi et al. (2017) in Section VII. Section VIII concludes the work.

II Description of a Molecular Hyperfine Interferometry Experiment

A surface-sensitive molecular hyperfine interferometer uses a beam of molecules to probe various surface properties. To do this, a set of magnetic fields are used to simultaneously manipulate the internal hyperfine states of the probe molecules and create a spatial superposition of molecular wavepackets. These wavepackets sequentially impact the sample surface and scatter in all directions. A second set of magnetic fields collects the molecules scattered in a narrow solid angle. This second set of magnetic fields further manipulates the molecular wavepackets, partially recombining them and allowing for molecular self-interference. Wavepackets with particular hyperfine states are then passed into a detector. A schematic of the experiment is depicted in Figure 1. We now discuss the different stages of the experiment in more detail.

The beam source must produce a continuous (or pulsed) beam of molecules with a sufficiently narrow velocity profile, mean velocity suitable for a particular experiment, sufficiently high flux, and a density low enough to ensure that the molecules are non-interacting. One current apparatus Godsi et al. (2017) uses a supersonic expansion to produce such a beam. One can also envision experiments with slow molecular beams produced by extraction (sometimes with hydrodynamic enhancement) from a buffer-gas cooled cell Hutzler et al. (2012) or with molecular beams controlled by electric-field Bethlem et al. (2002) or magnetic-field deceleration Wiederkehr et al. (2012). Deceleration provides control over the mean velocity and narrows the velocity spread Krems (2018), which could be exploited for novel interferometry-based applications.

The experiment selects molecules in particular hyperfine states by employing a magnetic lens whose magnetic field has a gradient in the radial direction. A cylindrically-symmetric field gradient is used to ensure sufficient molecular flux. The lens focusses molecules with low-field seeking states and defocusses molecules with high-field seeking states, allowing for purification of the molecular beam. After the lens, a unique quantization axis for the internal states is developed by using an auxiliary field that adiabatically rotates all magnetic moments until they lie along a single direction perpendicular to the beam path. The end of this auxiliary field is a strong dipolar field aligned along the $z$ direction that defines the quantization axis. Hexapole magnets can be used as a magnetic lens as their magnetic field gradients are sufficiently cylindrically symmetric Jardine et al. (2001); Dworski et al. (2004); Godsi et al. (2017). More details about the internal states of the molecules immediately after the magnetic lens can be found in Section III.

Solenoids whose magnetic fields are parallel to the beam propagation path are used to manipulate the molecular hyperfine states. These solenoids are helically-wrapped wire coils whose corresponding magnetic fields are generated by an electric current passing through each coil. These solenoids are labelled as the control fields in Figure 1. Arbitrary magnetic field profiles can be obtained by changing the solenoid winding patterns and/or using multiple successive solenoids.

The hyperfine states of a molecule change energy as the molecule enters a magnetic field. These changes to the hyperfine energy levels cause simultaneous changes in the molecular momenta, as the total energy is conserved. That is, when molecules enter a solenoid, molecules in low-field seeking states slow down and those in high-field seeking states speed up. Furthermore, because the direction of the magnetic field in a control field is not along the $z$ axis, the molecules are in a superposition of hyperfine states, with respect to the quantization axis defined by the magnetic field. Thus, the differences in momenta cause the different components of each molecular wavepacket to spatially separate as the wavepacket traverses the solenoid. Upon exiting the solenoid, the components of each wavepacket return to their original momenta, but remain spatially separated. That is, each wavepacket is now in an extended spatial superposition.

Each of these spatially separated wavepacket components comprise a superposition of the field-free hyperfine states. The exact superpositions of each wavepacket component, as well as the spatial separations between the components, depend on the magnetic field profile of the first branch. Each of the wavepacket components sequentially impacts the sample surface and scatters in all directions. However, the experiment only captures those molecules that pass through a particular solid angle. While a current experiment Godsi et al. (2017) fixes the angle between the two branches, one can in principle explore many different scattering geometries by varying both the angle between the two branches of the apparatus and the orientation of the sample.

After scattering, the collected molecules enter another set of control fields in the second branch of the apparatus. The hyperfine states again change in energy and momenta. In a helium-3 spin echo experiment, if the second magnetic field profile is identical but opposite in direction to the magnetic field profile of the first branch, the spatially separated wavepacket components realign (to first order) as they traverse the magnetic field(s), producing a spin echo. This allows the wavepacket components to interfere with each other. Interestingly, it has recently been shown Litvin et al. (2019) that echoes are also produced when the device operates with the fields in the same direction. With an arbitrary hyperfine Hamiltonian, such a realignment is only partial, though still useful. Experiments can be performed that explore either this spin-echo region or different relationships between the two magnetic field profiles, which may allow for a variety of insights about the sample surface. For example, the two field profiles can be different or the field magnitudes can be varied simultaneously, keeping $B_{1}=-B_{2}$. These different regimes of operation may produce different echoes, which can be collectively analyzed to provide more insight into molecule-surface interactions.

Additionally, as the spatially separated wavepacket components hit the surface sequentially, rather than simultaneously, any temporal changes in the surface that are on the time scale of the impact-time separation can differentially impact the phases of each wavepacket component. This may result in different interference patterns or even loss of coherence. This loss of coherence is the basis for the sensitivity of HeSE measurements to surface motion Alexandrowicz and Jardine (2007). Here, as in the recent experiment by Godsi et al. Godsi et al. (2017), we focus on surfaces whose dynamics are either much faster or much slower then the molecule-surface or wavepacket-surface interaction times. Note, however, that the current framework is suitable for extension to interaction regimes where the surface dynamics are comparable to these time scales.

After leaving the last solenoid of the second branch, the wavepackets pass through another auxiliary field that begins with a strong dipolar field in the $z^{\prime}$ direction. The auxiliary field then adiabatically connects magnetic moments aligned along the quantization axis to the radial direction of the final hexapole lens. This hexapole lens then focusses wavepackets with low-field seeking hyperfine states into the ionization detector and defocusses the rest. Finally, the ionization detector produces a current that is proportional to the molecular flux into the detector port. We describe how to calculate the molecular flux that enters the detector port in Section IV and the related transfer matrix formalism in Section V.

Analyzing the detector current as a function of the magnetic field profiles, the apparatus geometry, and the sample orientation can provide information about the interaction of the molecules with the sample surface. We discuss one possible analysis scheme in Section VI.

II.1 Molecular Hyperfine Hamiltonian

In principle, the only requirement for a molecular species to be suitable for molecular hyperfine interferometry is that the molecule have internal degrees of freedom whose energies are magnetic-field dependent. Such a requirement could be fulfilled by the presence of a nuclear spin, a rotational magnetic moment, or even an electronic spin. In practice, however, if the energy dependence on the magnetic field is too weak relative to the kinetic energy, state selection and state manipulation is difficult. On the other hand, if the dependence is too strong, the molecules may be difficult to control. Given these restrictions, we deem molecules that have a closed shell and are in an electronic state with zero orbital angular momentum to be most suitable for molecular hyperfine interferometry. In this case, the dominant interactions induced by magnetic fields are due to the nuclear magnetic spins of the molecules.

The hyperfine states of such a closed shell molecule with zero orbital angular momentum arise from coupling between the nuclear spin and the rotational degrees of freedom. Interactions of these hyperfine states with a magnetic field arise from the response of the nuclear and rotational magnetic moments to the external magnetic field. We assume that the hyperfine Hamiltonian, also referred to here as the Ramsey Hamiltonian Ramsey (1952), is of the following form:

[TABLE]

where $\vec{B}$ is the vector of the external magnetic field, assumed to be uniform across the molecule; $\hat{I}$ and $\hat{J}$ are the nuclear spin and rotational angular momentum operators, respectively; $I$ and $J$ are the nuclear spin and rotational angular momentum quantum numbers, respectively; $U$ contains all spin-rotational couplings (such as $\hat{I}\cdot\hat{J}$ or $\hat{I}^{2}\hat{J}^{2}$); $V$ contains all interactions of the nuclear spins with the magnetic field (such as $\hat{I}\cdot\vec{B}$); and $Q$ contains all interactions of the rotational angular momentum with the magnetic field (such as $\hat{J}\cdot\vec{B}$). Both $V$ and $Q$ are assumed to be proportional to positive powers of $|\vec{B}|$.

At large magnetic fields, $V$ and $Q$ dominate, making the eigenbasis $\ket{Im_{I}Jm_{J}}$, where $m_{I}$ and $m_{J}$ are the projections of the angular momenta $\vec{I}$ and $\vec{J}$ onto the external magnetic field direction, respectively. At zero field, $\hat{H}^{\text{R}}$ is diagonalized by $\ket{IJFM}$, where $\hat{F}=\hat{I}+\hat{J}$ is the total angular momentum operator and $M$ is the projection of $\vec{F}$ onto a chosen quantization axis. At intermediate fields, the eigenbasis is a function of the magnetic field and can be represented as a superposition of either $\ket{IJFM}$ or $\ket{Im_{I}Jm_{J}}$ states. Note that $M$ is a good quantum number at all field strengths. We call an eigenstate of $\hat{H}^{\text{R}}$ a Ramsey state, which we denote as $\ket{R}$ and which has the energy $E_{R}$. The number of eigenstates of $\hat{H}^{\text{R}}$ is $N_{R}$, such that $1\leq R\leq N_{R}$.

We treat the apparatus as a one dimensional system and account for the actual geometry by rotating the basis of the hyperfine states at the appropriate locations (see Section V.2). The total Hamiltonian can thus be written as

[TABLE]

where $\hat{p}$ is the centre of mass momentum operator, $m$ is the molecular mass, and $x$ is the position of the molecule in the apparatus. The magnetic field $\vec{B}(x)$ is now spatially dependent, reflecting the magnetic field profiles of the two branches of the apparatus.

In principle, the total Hamiltonian should incorporate molecule-surface interaction terms, such as the molecule-surface interaction potential. However, instead of treating the molecule-surface interactions explicitly, we include the interactions effectively through the use of a scattering transfer matrix (see Section V.3). This allows us to separate the details of the molecule-surface interaction from the propagation of the molecules through the apparatus. We can then treat the molecular propagation analytically while allowing for the scattering matrix to be determined by the level of theory practical for a particular system. Even more importantly, this approach allows us to treat the scattering matrix elements as free parameters that can be determined by fitting the calculated signal to an experimental signal. For the present manuscript, we treat the scattering matrix elements as arbitrary parameters, focussing primarily on the development of a theoretical formalism to describe the molecular propagation. We choose particular values for the scattering matrix elements only when we apply the formalism specifically to oH2 (Section VI). We also assume that the surface is static on timescales relevant to the experiment, such that the scattering matrix is time independent.

The system eigenstates $\ket{ER}$ are defined by the total Hamiltonian (2) through $\hat{H}\ket{ER}=E\ket{ER}$. Note that the system eigenstate $\ket{ER}$ is $N_{\text{R}}$ degenerate and that any linear combination of these states with the same label $E$ is also an eigenstate of $\hat{H}$. This degeneracy occurs as, while the $N_{\text{R}}$ different Ramsey states may have different energies, the kinetic energy can always be selected to maintain the same total energy. For the sake of convenience, we choose the orthonormal basis to be that defined by $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(x)\boldsymbol{\left.\right)}\ket{ER}=E_{R}\ket{ER}$ for $x\leq 0^{-}$. The zero of $x$ is defined to be immediately after the magnetic lens, while $y^{\pm}\equiv\lim_{\delta\to y^{\pm}}\delta$. We use these limit definitions as we will deal with discontinuities in the magnetic field when working with the transfer matrix formalism (Section V). As an example of the use of this notation, the statement that both one-sided limits are equal at the point $x$ $\left(i.e.\leavevmode\nobreak\ \lim_{a\to x^{-}}f(a)=\lim_{b\to x^{+}}f(b)\right)$ can be written as $f(x^{-})=f(x^{+})$.

The above definition of $\ket{ER}$ produces, for all $x$, a unique labelling of the system eigenstate $\ket{ER}$ by the total energy $E$ and the internal state $R$, where $R$ is a Ramsey state in the high magnetic field located immediately after the magnetic lens (i.e. at $x=0^{-}$). Note that, because of this definition, $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(x)\boldsymbol{\left.\right)}\ket{ER}\neq E_{R}\ket{ER}$ for $x\geq 0^{+}$; that is, the system eigenstates are superpositions of the local Ramsey states for $x\geq 0^{+}$. This unique labelling of the system eigenstates is valid for all $x$ as the eigenstate wavefunctions have a well-defined phase relationship throughout the entire apparatus. See Section V.1 for more details on the specifics of this phase relationship.

III Impact of the Magnetic Lens on the Molecular States

The magnetic lenses are designed to focus molecules with certain hyperfine states either onto the sample or into the detector. The remaining molecules are either defocussed or insufficiently focussed and contribute significantly less to the experimental signal. Roughly, high-field seeking states are defoccussed, some of the low-field seeking states are well focussed and the rest of the low-field seeking states are partially focussed. The actual proportions of each hyperfine state in the molecular beam must be measured or calculated from simulation. These magnetic lenses typically use large magnetic fields and large magnetic field gradients to perform this focussing Dworski et al. (2004); Jardine et al. (2001).

In general, magnetic lenses may take different forms, but we will consider lenses that have one key feature: the internal degrees of freedom of the outgoing molecular wavepackets are decohered in the high-magnetic field basis (i.e. $\ket{Im_{I}Jm_{J}}$). More precisely, we assume that the wavepacket exiting the magnetic lens is a mixed state of the form:

[TABLE]

where

[TABLE]

$\psi_{R_{0}k_{0}}(r)\equiv\braket{r}{\Psi_{R_{0}k_{0}}}$ is the wavefunction of a molecule in state $\ket{R_{0}}$; $\rho_{0}$ is the initial (time $t=0$) density matrix; $\ket{rR_{0}}\equiv\ket{r}\ket{R_{0}}$; $\ket{r}$ is an eigenstate of the position operator; $\ket{R_{0}}$ is an eigenstate of $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{lens}})$; $\vec{B}_{\mathrm{lens}}$ is a high magnitude, $z$-aligned magnetic field; $k_{0}$ is the experimentally-determined mean wavenumber of the wavepacket; and $P_{R_{0}}$ is the probability that the hyperfine statevector of the molecule is $\ket{R_{0}}$. Note that $\rho_{0}$ is diagonal in $\ket{R_{0}}$ but not in $\ket{r}$ (or $\ket{k}$, the momentum basis). Also, $\vec{B}_{\mathrm{lens}}=\vec{B}\left(x=0^{-}\right)$ corresponds to the final portion of the auxiliary field (i.e. a strong, $z$ aligned, dipolar field), not the field inside the hexapole magnet itself (see Section II).

That such a form of the wavepacket is valid follows from the work by Utz et al. Utz et al. (2015). The authors show that the two wavepackets arising from a spin–$\frac{1}{2}$ particle passing through a Stern-Gerlach apparatus are quickly decohered with respect to one another, even before they separate spatially. That is, the quantum dynamics themselves cause decoherence between the spin degrees of freedom (but not the spatial); a measurement or coupling to an external bath is not required. This decoherence occurs as the large magnetic field gradients cause a rapid oscillation in the off-diagonal terms of the extended Wigner distribution. That is, the phase relationship between the spin-up and spin-down components oscillates heavily in both the position and momentum bases, destroying coherence.

Given that the magnetic lenses we consider act like a Stern-Gerlach apparatus for the molecular hyperfine states, it is reasonable to assume that the internal hyperfine degrees of freedom will also decohere. Thus, we need only determine the values of $P_{R_{0}}$ for a specific magnetic lens. These can be found via semi-classical calculations Godsi et al. (2017); Krüger et al. (2018), may be measured experimentally Krüger et al. (2018) or may potentially be determined by solving the full 3D Schrödinger equation within the lens.

The mean velocity $v_{0}$ and velocity spread $\sigma_{v}$ of the molecules in the molecular beam can be measured experimentally Godsi et al. (2017). Both of these values are determined from the position and profile of scattering peaks obtained from the scattering of the probe molecules by appropriate sample surfaces Godsi et al. (2017). We assume that the initial wavefunction of a molecule $\psi_{R_{0}k_{0}}(r)$ is Gaussian and is characterized by $k_{0}\equiv mv_{0}/\hbar$ and $\sigma_{k}\equiv m\sigma_{v}/\hbar$, where $m$ is the mass of the molecule. More precisely,

[TABLE]

where $k_{0}^{R_{0}}$ is taken to be $k_{0}$. Though $k_{0}^{R_{0}}$ may in fact depend slightly (on the order of ppm) on $R_{0}$, we show later that the experimental signal is insensitive to small changes in $k_{0}^{R_{0}}$.

IV Wavepacket Propagation and Signal Calculation

The primary measured value of the experiment is a current that is proportional to the molecular flux entering the detector. This measured current is a function of the magnetic fields, the scattering geometry, and the surface properties. The molecular flux entering the detector can be calculated as the product of the molecular flux incident to the apparatus and the probability that a molecule entering the apparatus will successfully pass through the apparatus and be detected. It is this probability of detection $P_{\text{detection}}$ that is sensitive to the experimental parameters and surface properties. Note that the incident molecular flux could be either continuous or pulsed, as long as the density is low enough that the molecules can be considered non-interacting.

As the detector has a finite time-response, the probability of detection is given by

[TABLE]

where $t_{1}$ and $t_{2}$ are the initial and final times of the detection window $\tau=t_{2}-t_{1}$, and $\langle\hat{C}(t)\rangle$ is the expectation value of the detector measurement operator $\hat{C}$. This expectation value is given by

[TABLE]

where $\hat{\rho}(t)\equiv\hat{U}\rho_{0}\hat{U}^{\dagger}=\sum_{R_{0}}P_{R_{0}}\ket{\Psi_{R_{0}k_{0}}(t)}\bra{\Psi_{R_{0}k_{0}}(t)}$ is the time evolved density matrix, $\hat{U}\equiv e^{-i\frac{\hat{H}}{\hbar}t}$ is the time evolution operator, $\rho_{0}$ is the density matrix (3) at $t=0$, and $\ket{\Psi_{R_{0}k_{0}}(t)}\equiv\hat{U}\ket{\Psi_{R_{0}k_{0}}}$.

Given that the detector consists of a magnetic lens that focusses molecules with particular states into a measuring apparatus, such as an ionization detector Godsi et al. (2017), and that the internal degrees of freedom of these molecules are decohered by the second magnetic lens (see Section III), we can model the detector with a diagonal operator

[TABLE]

The matrix elements of $\hat{C}$ are the probabilities $c_{R_{D}}(x)$ of detecting, at position $x$, a molecule whose internal state is a high-field eigenstate $\ket{R_{D}}$ of $\hat{H}^{\text{R}}$. Note that $c_{R_{D}}(x)=0\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ x<x_{D}$, the detector position.

Using the time evolution operator, we determine the time-dependence of the density matrix $\rho(t)$ to be

[TABLE]

where $\alpha_{k_{0}R_{0}}^{ER}\equiv\int\text{d}r\leavevmode\nobreak\ \psi_{R_{0}k_{0}}(r)\Phi_{R_{0}}^{*ER}(r)$ is the overlap between the initial wavefunction $\psi_{R_{0}k_{0}}(r)$ and the system eigenstate wavefunction $\Phi_{R_{0}}^{ER}(r)\equiv\braket{rR_{0}}{ER}$.

We can evaluate $\langle\hat{C}(t)\rangle$ by inserting a resolution of the identity $\sum_{R_{D}}\int\text{d}r\leavevmode\nobreak\ \ket{rR_{D}}\bra{rR_{D}}$, where $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(x)\boldsymbol{\left.\right)}\ket{R_{D}}=E_{R_{D}}\ket{R_{D}}$ for $x\geq x_{D}^{+}$ and $x_{D}$ is the starting location of the detector (see Figure 2). In other words, $\ket{R_{D}}$ is a Ramsey state in the strong dipolar magnetic field of the detector auxiliary field. The result is

[TABLE]

where we have evaluated the trace in the $\ket{r^{\prime}R_{D}^{\prime}}$ basis and

[TABLE]

We emphasize that $R$ and $R_{D}$ are indices of different sets of Ramsey states, i.e. $\braket{R}{R_{D}}\neq\delta_{RR_{D}}$, unless the magnetic fields at the first magnetic lens ($x=0^{-}$) and the detector magnetic lens ($x=x_{D}^{+}$) happen to be identical.

We also have

[TABLE]

which, when inserted with Eqn. (11) into Eqn. (10), results in

[TABLE]

The initial wavepacket is almost entirely confined to the region $r\leq 0^{-}$, as $\psi_{R_{0}k_{0}}(r)$ has a Gaussian profile (5) with spatial width on the order of $10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (as determined from the measured velocity distribution for oH2 Godsi et al. (2017)). Thus, we can evaluate $\alpha_{k_{0}R_{0}}^{ER}\equiv\int\text{d}r\leavevmode\nobreak\ \psi_{R_{0}k_{0}}(r)\Phi_{R_{0}}^{*ER}(r)$ if we know $\Phi_{R_{0}}^{ER}(r)$ for $r\leq 0^{-}$. Given the definition of the eigenstate $\ket{ER}$, discussed in Section II.1, we show in Section V.1 that $\Phi_{R_{0}}^{ER}(r)=A_{R}e^{irk^{ER}}\delta_{RR_{0}}$ for $r\leq 0^{-}$ (cf. Eqn. (34)), where $k^{ER}\equiv\frac{\sqrt{2m\left(E-E_{\text{R}}\right)}}{\hbar}$ (cf. Eqn. (23)). Combined with the definition (4) of $\psi_{R_{0}k_{0}}(r)$,

[TABLE]

where $\Gamma^{ER}_{k_{0}R_{0}}=A_{R}^{*}\frac{(2\pi)^{\frac{1}{4}}}{\sqrt{\sigma_{k}}}e^{-\frac{\left(k^{ER}-k_{0}^{R_{0}}\right)^{2}}{4\sigma_{k}^{2}}}$. Thus,

[TABLE]

where we have performed the sums over $R$ and $R^{\prime}$.

If the detection window $t_{2}-t_{1}$ is large enough that the entire wavepacket passes through the detection region defined by $c_{R_{D}}(z)$ we have

[TABLE]

where $\frac{2\pi\hbar}{\tau}\delta\left(E-E^{\prime}\right)=\frac{1}{\tau}\int_{-\infty}^{\infty}\text{d}te^{-\frac{i}{\hbar}(E-E^{\prime})t}$.

Physically, one can see that the probability of detection is proportional to the overlap $\left|\Gamma^{ER_{0}}_{k_{0}R_{0}}\right|^{2}$ of the initial wavepacket and a system eigenstate multiplied by the overlap $\int\text{d}r\leavevmode\nobreak\ \frac{2\pi\hbar}{\tau}\left|\Phi_{R_{D}}^{ER_{0}}(r)\right|^{2}c_{R_{D}}(r)$ of the same system eigenstate and the detection region, as expected.

Substituting for $\left|\Gamma^{ER_{0}}_{k_{0}R_{0}}\right|^{2}$ and given that

[TABLE]

for $r\geq x_{D}^{+}$ (cf. Eqn. (34)), we have

[TABLE]

where $c_{R_{D}}\equiv\int\text{d}r\leavevmode\nobreak\ c_{R_{D}}(r)$ and $\beta_{R_{D}}^{ER_{0}}\equiv\braket{R_{D}}{ER_{0}}$, the projection of the system eigenstate $\ket{ER_{0}}$ onto the detector eigenstate $\ket{R_{D}}$ at $x_{D}^{+}$. For the purposes of comparing to experiment, only the dependence of $P_{\text{detection}}$ on the experimental parameters is needed, not its absolute value. Also, the value of $A_{R_{0}}=1$ as $A_{R}e^{irk^{ER}}\equiv\braket{rR_{0}}{ER_{0}}=e^{irk^{ER}}(\mathrm{for}\leavevmode\nobreak\ r\leq 0^{-})$ because of the specific definition of the system eigenstates (see Section II.1). Additionally, one can see that $P_{\text{detection}}$ is not sensitive to minor (on the order of ppm) changes in $k_{0}^{R_{0}}$ as $\sigma_{k}\propto k_{0}$ in experiment Godsi et al. (2017). Finally, in Eqn. (17), only $\beta_{R_{D}}^{ER_{0}}$ is dependent on the magnetic fields, the scattering geometry, and the surface properties. It is thus sufficient to work with the following equation:

[TABLE]

To determine the values of $\beta_{R_{D}}^{ER_{0}}$, we derive and apply the transfer matrix method with internal degrees of freedom (Section V).

V Transfer Matrix Formalism with Internal Degrees of Freedom

The transfer matrix method as applied in quantum transport turns the solution of the time-independent Schrödinger equation of a 1D system into a product of matrices Walker and Gathright (1992). Pedagogical introductions can be found in Refs. Walker and Gathright (1992); Sánchez-Soto et al. (2012); Mello and Kumar (2004). The present problem has two unique features: (i) the propagating molecules have many internal degrees of freedom which may be mixed as the molecule transitions from one local field to another and (ii) molecules change their propagation direction after scattering by the surface. Problem (i) is addressed in Section V.1, while (ii) is addressed in Section V.2. The impact of scattering on the internal degrees of freedom is accounted for by using a scattering transfer matrix (Section V.3).

The transfer matrix formalism we present in Section V.1 is similar to the mixed multicomponent transfer matrix formalism described in Diago-Cisneros et al. (2006) and can be viewed as an extension and application of the transfer matrix formalism used in the study of molecular tunnelling Saito and Kayanuma (1994); Jarvis and Bulte (1998). The formalism combines transfer matrices that incorporate the internal molecular degrees of freedom of a composite particle Saito and Kayanuma (1994); Jarvis and Bulte (1998) with eigenbasis changes between regions of the external potential. Similar eigenbasis changes have been employed in the transfer matrix formalism used in the envelope function approximation, which is used to calculate electronic properties in abrupt semi-conductor heterostructures Foreman (1996). We further extend the transfer matrix formalism in Sections V.2 and V.3 to account for the impact of scattering on the molecules and their relevant internal degrees of freedom.

V.1 Propagation and Discontinuity Matrices

We first break up the arbitrary magnetic field profiles of the apparatus into rectangular regions of constant field, as shown in Figure 2. We then solve the Schrödinger equation for a single eigenstate in a single region of constant field. Subsequently, we determine the impact of the boundary conditions that exist at the discontinuity between two regions of constant field. Using these solutions, we determine matrices that describe the spatial dependence of the eigenstate wavefunction coefficients within a region of constant field (propagation matrices) and matrices that describe how these coefficients change across the discontinuity between two regions of constant field (discontinuity matrices). Note that while we derive these matrices for molecules whose internal degrees of freedom are described by the Ramsey Hamiltonian (1), the formalism is not limited to this Hamiltonian.

Within a region of uniform magnetic field, the Ramsey Hamiltonian $\hat{H}^{\text{R}}$ is constant, which allows us to derive the propagation matrix that includes the internal degrees of freedom. We begin by expanding a system eigenstate $\ket{E\tilde{R}}$ as

[TABLE]

where $\Phi_{R}^{E\tilde{R}}(x)\equiv\braket{xR}{E\tilde{R}}$, we define $\ket{xR}\equiv\ket{x}\ket{R}$, and $\ket{R}$ is one of the $N_{R}$ Ramsey states of a molecule in some magnetic field $\vec{B}$. Note that $\vec{B}$ is not necessarily the local magnetic field $\vec{B}_{\mathrm{loc}}$ of the current region and thus $\ket{R}$ is not necessarily an eigenstate of $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})$ at this point. Also, the eigenstates $\ket{E\tilde{R}}$ are labelled by their energy $E$ and a particular Ramsey index $\tilde{R}$, such that \hat{H}^{\text{R}}(\underaccent{\tilde}{\vec B})\ket{E\tilde{R}}=E_{\tilde{R}}\ket{E\tilde{R}}, with \underaccent{\tilde}{\vec B} an arbitrarily chosen magnetic field.

Using Eqn. (19), the Schrödinger equation with the total Hamiltonian (2) can be shown to be (Appendix A):

[TABLE]

where $H^{\text{R}}_{R_{0}R}=\bra{R_{0}}\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})\ket{R}$. Eqn. (20) is in general difficult to solve because of the coupling of the internal degrees of freedom by $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})$. However, if we choose the eigenbasis of the internal degrees of freedom to satisfy $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})\ket{R}=E_{\text{R}}\ket{R}$ (that is, $\ket{R}$ is now a Ramsey state of a molecule in the local magnetic field $\vec{B}_{\mathrm{loc}}$), the equations decouple and we obtain

[TABLE]

The solution is

[TABLE]

where $A_{R}$ and $B_{R}$ are $R$-dependent coefficients and

[TABLE]

As per the single channel transfer matrix method Walker and Gathright (1992), given that $\Phi_{R}^{E\tilde{R}}\left(x+\Delta x\right)=A_{R}e^{ik_{R}x}e^{ik_{R}\Delta x}+B_{R}e^{-ik_{R}x}e^{-ik_{R}\Delta x}$, we can collect the $A_{R}$ and $B_{R}$ coefficients into a $2N_{R}$-dimensional coefficient vector $\vec{\phi}_{x}=\left(A_{1},A_{2},...,A_{N_{R},}B_{1},B_{2},...,B_{N_{R}}\right)^{T}$ and write

[TABLE]

where $\mathbf{\Pi}_{x}$ is the $2N_{R}\times 2N_{R}$ propagation matrix

[TABLE]

where $\oplus$ denotes the direct sum.

Following the derivation of Ref. Walker and Gathright (1992), we can determine how the coefficients transform across a step discontinuity in the magnetic field. Using the propagation matrix (25) and a relabelling of the coordinate system, we can always set the discontinuity to appear at $x=0$. Given that Eqn. (20) applies everywhere, the coefficients $\Phi_{R}^{E\tilde{R}}(x)$ and their derivatives are continuous across the discontinuity (i.e. $\Phi_{R}^{E\tilde{R}}(x)\in C^{1}(x)$), for each value of $R$. However, the coefficients are only known when $\ket{R}$ is an eigenstate of $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})$, which differs on each side of the discontinuity (that is, $\vec{B}(0^{-})\neq\vec{B}(0^{+})$). Note that the wavevector $\ket{E\tilde{R}}$ is the same everywhere in the system. Thus, by writing the wavevector $\ket{E\tilde{R}}$ in the two different bases corresponding to the eigenstates of $\hat{H}^{\text{R}}$ on each side of the field, we see that the coefficients at a specific value of $x$ are related by a basis transformation:

[TABLE]

where $\ket{E\tilde{R}^{\pm}}$ is the wavevector written in the basis of $\ket{R^{\pm}}$, $\ket{R^{\pm}}$ are the eigenstates of $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{\pm})\boldsymbol{\left.\right)}$ on the left $(-)$ and right $(+)$ sides of the discontinuity at $x=0$, respectively, and $\sum_{R^{+}}\ket{R^{+}}\bra{R^{+}}$ was inserted in the third line (recall that $\ket{xR^{-}}\equiv\ket{x}\ket{R^{-}}$). The values $\braket{R^{-}}{R^{+}}$ are recognized as the matrix elements $S_{R^{-}R^{+}}$ of the matrix $\mathbf{S}^{R^{+}}_{R^{-}}$ whose columns are the eigenstates of $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{+})\boldsymbol{\left.\right)}$ written in the $\ket{R^{-}}$ basis.

Since $\Phi_{R}^{E\tilde{R}}(x)\in C^{1}(x)$ for each value of $R$ separately, we can equate the two limits $\lim_{x\to 0^{\mp}}\Phi_{R^{+}}^{E\tilde{R}}(x)$ and the two limits of the derivative $\lim_{x\to 0^{\mp}}\frac{\partial}{\partial x}\Phi_{R^{+}}^{E\tilde{R}}(x)$. Solving the resultant equations for the coefficients $A_{R^{+}}$ and $B_{R^{+}}$, we obtain (Appendix B):

[TABLE]

where $S_{R^{-}R^{+}}^{*}\equiv\braket{R^{+}}{R^{-}}$, $\Delta_{R^{+}R^{-}}^{\pm}\equiv\frac{1}{2}\left(1\pm\frac{k_{R^{-}}}{k_{R^{+}}}\right)$, $k_{R^{\pm}}\equiv\frac{\sqrt{2m\left(E-E_{\text{R}^{\pm}}\right)}}{\hbar}$, and $E_{\text{R}^{\pm}}\equiv\braket{R^{\pm}}{\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{\pm})\boldsymbol{\left.\right)}}{R^{\pm}}$. There are $N_{R}$ such sets of equations, one for each value of $R^{+}$. Working again with $\vec{\phi}_{x}=\left(A_{1},A_{2},...,A_{N_{R},}B_{1},B_{2},...,B_{N_{R}}\right)^{T}$, one can write the matrix equation

[TABLE]

where $x^{\mp}$ indicates the location just before $({-})$ or just after $({+})$ the discontinuity located at $x$ and $\mathbf{K}$ is the $2N_{R}\times 2N_{R}$ discontinuity matrix

[TABLE]

where $\circ$ denotes the element-wise Hadamard product, such that $({\mathbf{S}^{R^{+}}_{R^{-}}}^{\dagger}\circ\mathbf{\Delta^{\pm}})_{R^{+}R^{-}}\equiv S_{R^{-}R^{+}}^{*}\Delta_{R^{+}R^{-}}^{\pm}$. This matrix allows one to calculate the coefficients of the wavefunction as one moves from one region of constant magnetic field to another through a discontinuity. Thus, if one breaks up any magnetic field profile into a series of constant regions separated by discontinuities, one can systematically approach a perfect description of the propagation of a molecule with internal degrees of freedom through a magnetic field of arbitrary profile through repeated application of $\mathbf{K}$ and $\mathbf{\Pi}_{x}$. Furthermore, this approach is not restricted to molecules moving through magnetic fields. Many other types of quantum objects moving in a single dimension with internal degrees of freedom that couple to an external static potential can also be analyzed in this way.

The above analysis indicates that one needs to keep track of $2N_{R}$ components to build up the eigenstates of the system exactly. However, for the current application in mind, one only needs $N_{R}$ components as the magnetic fields typically change the linear molecular momentum by such a small amount that the amplitudes $B_{R}$ of the reflected part of the wavefunction are negligible. That is, any backscattering of the molecules by the magnetic fields is negligible and can be ignored.

For example, a typical velocity of the oH2 molecules in the experiment of Ref. Godsi et al. (2017) is $v_{H_{2}}=$1450\text{,}\mathrm{m}\mathrm{/}\mathrm{s}$$. This corresponds to the kinetic energy $E_{H_{2}}=\frac{1}{2}m_{H_{2}}v_{H_{2}}^{2}=5.31\times 10^{9}\text{ kHz}$. The data reported by Ramsey Ramsey (1952) indicates that the maximum energy change for the hyperfine states of oH2 at 500G is approximately -2550 kHz. The experiment of Ref. Godsi et al. (2017) has magnetic fields up to about 1000G. For such fields, the energy changes are approximately linear, so we expect the maximum change in energy to be $\Delta E\approx-5100$ kHz. In the field-free region before the discontinuity, $k_{R^{-}}\approx m_{H_{2}}v_{H_{2}}/\hbar$ and after the discontinuity in the field, $k_{R^{+}}\approx\sqrt{2m_{H_{2}}(E-\Delta E)}/\hbar$, as per Eqn. (23). Then, $|\Delta_{R^{+}R^{-}}^{-}|\approx 2.4\times 10^{-7}$ and $|\Delta_{R^{+}R^{-}}^{+}|\approx 1$, making $\mathbf{K}$ approximately diagonal and illustrating the decoupling of the forward and backward channels under typical experimental conditions.

We thus only need to keep track of the $A_{R}$ components, which correspond to the forward-propagating momenta. We can define a new coefficient vector

[TABLE]

The corresponding $N_{R}\times N_{R}$ propagation $\mathbf{P}_{x}$ and discontinuity $\mathbf{D}$ matrices are

[TABLE]

where the matrix elements of ${\mathbf{S}^{R^{+}}_{R^{-}}}^{\dagger}$ are $S_{R^{-}R^{+}}^{*}\equiv\braket{R^{+}}{R^{-}}$, $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{\pm})\boldsymbol{\left.\right)}\ket{R^{\pm}}=E_{R^{\pm}}\ket{R^{\pm}}$, $0^{\pm}$ indicates the position just to the left $({-})$ or right $({+})$ of the discontinuity, and $k_{R}$ is defined as in Eqn. (23). Specifically, $\mathbf{D}$ changes the basis of the vector $\vec{\psi}_{x}$ from $\ket{R^{-}}$ to $\ket{R^{+}}$. That is, $\vec{\psi}_{x}$ is always in the eigenbasis of $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(x)\boldsymbol{\left.\right)}$. Finally, given that $B_{R}\approx 0$, the eigenstate coefficients are now

[TABLE]

Given that a generic transfer matrix $\mathbf{M}$ has the property $\mathbf{M}\mathbf{\sigma}_{z}\mathbf{M}^{\dagger}=\mathbf{\sigma}_{z}$Walker and Gathright (1992), the decoupling of the forward and backward channels implies that the forward channel matrix $\mathbf{M}_{F}$ (composed of a product of $\mathbf{P}_{x}$ and $\mathbf{D}$ matrices) is now unitary.

V.2 Rotation Matrices

Scattering by the sample surface changes both the propagation direction and the internal states of the molecule. To take into account the change in the direction of the propagation path when applying the transfer matrix formalism, we need only change the orientation of the quantization axis. However, to address the impact of scattering on the internal states, we need to apply a scattering matrix that is written with respect to a particular reference frame (which is often a sample-fixed frame, see Section V.3). Thus, instead of just rotating the quantization axis from the first branch to the second branch (to account for the change in the direction of propagation), we need to first rotate from the initial reference frame ($xyz$ in Figure 1) to the reference frame of the scattering matrix. Then, after applying the scattering matrix, we need to rotate from the scattering matrix reference frame to the final reference frame ($x^{\prime}y^{\prime}z^{\prime}$ in Figure 1). To perform these rotations coherently, we apply $N_{\mathrm{R}}\times N_{\mathrm{R}}$ rotation matrices $\mathbf{R}(\phi,\Theta,\chi)$ to $\vec{\psi}_{x}$, where $\phi$, $\Theta$, and $\chi$ are the Euler angles in the $ZYZ$ convention (with $Y$ and $Z$ being the axes of a space-fixed frame; see Ref. Zare (1988)). In this way, we can account for both specular and non-specular scattering geometries and for various orientations of the sample surface.

To change the orientation of the quantization axis, we perform passive rotations on the state vector $\vec{\psi}_{x}$. These passive rotations modify the basis of $\vec{\psi}_{x}$, but leave the physical state unchanged. For example, if we were to assume that the only impact of scattering was to change the propagation direction, we would need to perform a passive rotation of the state vector about the $y$ axis by the angle $\theta$ to account for a change of angle $\theta$ in the propagation direction (for the definition of the axes shown in Figure 1). We would perform this rotation by applying the equivalent active rotation of angle $-\theta$ to $\vec{\psi}_{x}$; that is, by using the matrix $\mathbf{R}(0,-\theta,0)$.

For the general case, we work with the rotation matrices $\mathbf{R}(\phi,\Theta,\chi)$, whose matrix elements, when written in the $\ket{R}$ eigenbasis of $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})$ where $\vec{B}_{\mathrm{loc}}$ is the local magnetic field, are:

[TABLE]

where $\ket{FM}\equiv\ket{IJFM}$ is an angular momentum state with total angular momentum $F$, $z$ axis projection $M$, total nuclear spin angular momentum $I$ and total rotational angular momentum $J$; the subscripts of $\mathbf{R}_{R}$ and $\mathbf{R}_{FM}$ denote the basis of the matrix representation, $\ket{R}$ and $\ket{FM}$, respectively; $\mathbf{S}^{R}_{FM}$ is the matrix whose columns are the eigenstates $\ket{R}$ written in the $\ket{FM}$ basis, $\hat{R}(\phi,\Theta,\chi)$ is the rotation operator (with the same $ZYZ$ convention mentioned above) and

[TABLE]

where $D^{F}_{M^{\prime}M}(\phi,\Theta,\chi)$ are the Wigner D-matrices and $d^{F}_{M^{\prime}M}(\Theta)$ are the Wigner small d-matrices Zare (1988). Note that $\mathbf{R}_{FM}(\phi,\Theta,\chi)$ is diagonal in $F$, because of conservation of angular momentum, but not diagonal in $M$ Zare (1988). Thus, one must be careful to also perform a passive rotation on the local magnetic field vector if rotations are performed in a region with non-zero field. Typically, however, the sample chamber is magnetically shielded.

We also note that the rotation may impact how to appropriately match the boundary conditions between the eigenstate immediately after rotation and the eigenstate at the start of the second branch. As the propagation matrices $\mathbf{P}_{x}$ (32) are defined with respect to the momentum, which may be positive or negative, it is important to choose the sign of the momentum that results in the probability current flowing in the same direction as the molecular propagation. For example, using the axis definitions in Figure 1, $+k_{R}$ is chosen for the first branch and $-k_{R}$ for the second branch.

V.3 Scattering Transfer Matrices

Scattering by the sample surface can involve many complex phenomena that may change the internal state, the momentum, and the total energy of the scattering molecule. For the present manuscript, we focus on scattering processes that conserve the total energy of the molecules. Energy-conserving scattering may, however, include transfer of energy between the internal and translational degrees of freedom. Such scattering processess are described by a general, non-diagonal scattering matrix in the basis of the molecular states.

The interactions of the molecules with the sample surface can be phenomenologically described with the total scattering transfer matrix. This matrix is the $2N_{\mathrm{R}}\times 2N_{\mathrm{R}}$ matrix ${\mathbf{\tilde{\Sigma}}}$ that relates the wavefunctions on the “left” side of the scattering event to those on the “right” (as opposed to the scattering matrix, which relates the incoming wavefunctions to the outgoing). However, because the initial wavepacket (5) does not contain any negative momentum states, the magnetic fields of the solenoids do not cause significant backscattering (Section V.1), and the detector only detects molecular flux in the forward scattering direction, we need only work with the $N_{\mathrm{R}}\times N_{\mathrm{R}}$ matrix $\mathbf{\Sigma}\equiv\mathbf{P}_{\mathrm{fwd}}\mathbf{\tilde{\Sigma}}{\mathbf{P}^{\dagger}_{\mathrm{fwd}}}$, where $\mathbf{P}_{\mathrm{fwd}}$ is an $N_{\mathrm{R}}\times 2N_{\mathrm{R}}$ projection matrix onto the forward scattering states. We define $\mathbf{\Sigma}$ in the $\ket{Im_{I}Jm_{J}}$ basis, where the $\ket{Im_{I}Jm_{J}}$ states are themselves defined with respect to the quantization axis that is normal to the surface sample. We choose this basis to relate to scattering calculations, which are frequently carried out in the $\ket{Jm_{J}}$ basis with a quantization axis normal to the sample surface. In principle, however, any suitable set of Ramsey states $\ket{R_{\Sigma}}$ could be chosen as the basis for the scattering transfer matrix and any suitable quantization axis could be chosen, to take advantage of relevant symmetries.

In general, the scattering transfer matrix elements $\Sigma_{Im_{I}Jm_{J}I^{\prime}m^{\prime}_{I}J^{\prime}m^{\prime}_{J}}$ are functions of the incident energy $E$, the outgoing energy $E^{\prime}$, the incident momentum $\vec{k}$, and the outgoing momentum $\vec{k}^{\prime}$. As we are restricting ourselves to iso-energetic processes, $E=E^{\prime}$. Also, Eqn. (23) defines the magnitudes of the momentum before and after the scattering event. This leaves the scattering transfer matrix elements as functions of only energy and the four angles that define the scattering geometry. These angles are: (1) the angle between the two branches, (2) the angle between the surface normal and the scattering plane, (3) the angle between the first branch and the projection of the surface normal on the scattering plane, and (4) the azimuthal angle of the sample. The scattering plane is the plane defined by the two branches of the apparatus.

Given that the experiment only probes a single scattering direction at a time (see Section II and Figure 1), the scattering transfer matrix will not, in general, be unitary. This incorporates state-dependent loss channels into the formalism. Additionally, the scattering transfer matrix is, in general, time-dependent. Here, we assume that the time-scales of the surface dynamics are significantly different from the molecule-surface or wavepacket-surface interaction time-scales and assume $\mathbf{\Sigma}$ to be time-independent.

Because $\mathbf{\Sigma}$ is defined with respect to the surface normal, we use rotation matrices to appropriately change the basis of $\vec{\psi}$ before and after applying the scattering transfer matrix. We ensure that the total rotation corresponds to the change in propagation direction induced by scattering off of the sample surface and that the quantization axis is again coplanar with the two branches of the apparatus.

The scattering transfer matrix elements for a specific molecule-surface interaction can be determined from scattering calculations Godsi et al. (2017); Díaz et al. (2009). Alternatively, they can be treated as free parameters and determined from the experimental measurements by solving the inverse scattering problem. Such a problem can potentially be solved efficiently using machine learning based on Bayesian optimization Vargas-Hernández et al. (2019); Krems (2019).

V.4 Calculation of Eigenstate Coefficients

To determine the dependence of the probability of detection (18) on the magnetic fields and the surface properties, we must determine the coefficients $\beta_{R_{D}}^{ER_{0}}$. This can be done by multiplying the initial coefficient vector $\vec{\psi}^{ER_{0}}_{x_{0}}$ (31) of a system eigenstate $\ket{ER_{0}}$ by a succession of transfer matrices to obtain the final coefficient vector $\vec{\psi}^{ER_{0}}_{x_{D}}\equiv\left(\beta_{1}^{ER_{0}},\beta_{2}^{ER_{0}},\cdots,\beta_{N_{\mathrm{R}}}^{ER_{0}}\right)^{T}$:

[TABLE]

where $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ describe the propagation through the first and second branches of the apparatus, respectively, $\mathbf{M}_{\Sigma}$ describes the scattering, and ${\mathbf{S}^{R_{D}}_{R_{N}}}^{\dagger}$ changes the basis of the coefficient vector to the eigenbasis $\ket{R_{D}}$ of $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(x_{D}^{+})\boldsymbol{\left.\right)}$ at the location of the detector $x_{D}$. The $\mathbf{M}$ matrices are defined as

[TABLE]

where $R_{ini}$ refers to the eigenbasis $\ket{R_{ini}}$ of $\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{-})\boldsymbol{\left.\right)}$ at the initial location of the wavepacket; $R_{i}$ refers to the eigenbasis $\ket{R_{i}}$ of $\hat{H}^{\text{R}}(\vec{B}_{i})$ in region $i$ of the apparatus, as depicted in Figure 2; $FM$ refers to the $\ket{IJFM}$ basis where $\vec{F}\equiv\vec{I}+\vec{J}$ and $M$ is the projection on the local $z$ axis; $\alpha$, $\beta$, and $\gamma$ are the Euler angles that rotate the reference frame of the first branch ($xyz$ in Figure 1) onto the reference frame of the scattering transfer matrix, whose quantization axis is normal to the sample surface (see Section V.2 and Section V.3); $\alpha^{\prime}$, $\beta^{\prime}$, and $\gamma^{\prime}$ are the Euler angles that rotate the scattering transfer matrix reference frame onto the reference frame of the second branch ($x^{\prime}y^{\prime}z^{\prime}$ in Figure 1); $L_{i}$ is the signed length of region $i$, as depicted in Figure 2; the sign of $L_{i}$ indicates the direction of propagation with respect to the local $x$ or $x^{\prime}$ axis; $N$ is the total number of regions between $x_{0}$ and $x_{D}$ (see Figure 2); $n$ is the number of regions between the initial position of the wavepacket $x_{0}=0$ and the sample position $x_{S}$; $\mathbf{\Sigma}_{FM}\equiv{\mathbf{S}^{FM}_{R_{IJ}}}^{\dagger}\mathbf{\Sigma}{\mathbf{S}^{R_{IJ}}_{FM}}^{\dagger}$ is the scattering transfer matrix written in the $\ket{IJFM}$ basis; $R_{IJ}\equiv Im_{I}Jm_{J}$; and $\mathbf{\Sigma}$ is the scattering transfer matrix in the $\ket{Im_{I}Jm_{J}}$ basis. Note the product of the two rotation matrices $\mathbf{R}_{FM}(\alpha^{\prime},\beta^{\prime},\gamma^{\prime})\cdot\mathbf{R}_{FM}(\alpha,\beta,\gamma)=\mathbf{R}_{FM}(\phi,\Theta,\chi)$, where $\phi$, $\Theta$, and $\chi$ are the Euler angles that rotate the reference frame $xyz$ onto the frame $x^{\prime}y^{\prime}z^{\prime}$ (see Figure 1). All of the Euler angles mentioned above are in the $ZYZ$ convention, with $Y$ and $Z$ being the axes of a space-fixed frame and as per the convention defined in Ref. Zare (1988). Note also that while the scattering transfer matrix $\mathbf{\Sigma}$ is written here in the $\ket{Im_{I}Jm_{J}}$ basis, other suitable bases $\ket{R_{\Sigma}}$ may be used (see Section V.3), where $R_{\Sigma}$ refers to an arbitrary set of Ramsey states. In such a case, $\mathbf{\Sigma}_{FM}\equiv{\mathbf{S}^{FM}_{R_{\Sigma}}}^{\dagger}\mathbf{\Sigma}{\mathbf{S}^{R_{\Sigma}}_{FM}}^{\dagger}$. Also, note that the propagation matrices $\mathbf{P}_{L_{i}}$ (32) are defined with momentum $+k_{R}$ if the molecular propagation is in the direction of the local $x$ or $x^{\prime}$ axis or, conversely, with the momentum $-k_{R}$ if the molecular propagation is in the opposite direction of the local $x$ or $x^{\prime}$ axis (see Section V.2).

By defining a matrix $\mathbf{\Psi}^{E}_{x_{i}}\equiv\left(\vec{\psi}^{E1}_{x_{i}},\vec{\psi}^{E2}_{x_{i}},\cdots,\vec{\psi}^{EN_{\mathrm{R}}}_{x_{i}}\right)$, all $N_{\mathrm{R}}\times N_{\mathrm{R}}$ coefficients $\beta_{R_{D}}^{ER_{0}}$ can be simultaneously obtained from

[TABLE]

where $\mathbf{\Psi}^{E}_{x_{0}}\equiv\mathbb{1}_{N_{\mathrm{R}}}$ because of the specific definition of the system eigenstates (see Section II.1). Using Eqns. (38–41), we can obtain $\beta_{R_{D}}^{ER_{0}}$, and thus $P_{\text{detection}}$ (18), as functions of the magnetic field profile, the scattering matrix elements, and the scattering geometry.

VI Application to ortho-Hydrogen

The theoretical framework described in Sections II through V connects the scattering transfer matrix elements $\Sigma_{Im_{I}Jm_{J}I^{\prime}m^{\prime}_{I}J^{\prime}m^{\prime}_{J}}$ to the experimentally observed signal, which is proportional to $P_{\text{detection}}$ (18). By changing the magnetic field profiles in the two arms of the apparatus, one can obtain information about how the scattering affects various hyperfine states. To illustrate our theoretical framework and to demonstrate the impact of the scattering transfer matrix on the experimentally observed signal, we consider a beam of rotationally cold oH2 and a simplified apparatus that contains only a few regions of constant magnetic field, as depicted in Figure 3.

VI.1 Rotationally Cold Ortho-Hydrogen Hyperfine Hamiltonian

The Hamiltonian describing the relevant internal degrees of freedom of rotationally cold oH2 is Ramsey (1952):

[TABLE]

where, for simplicity, we have neglected magnetic shielding of the nuclear and rotational magnetic moments by the molecule and diamagnetic interactions of the molecule with the magnetic field; $\vec{B}$ is the local magnetic field; $\hat{I}$ is the nuclear spin operator; $\hat{J}$ is the rotational angular momentum operator; $\alpha\equiv\frac{\mu_{I}}{hI}\approx$4.258\text{,}\mathrm{k}\mathrm{H}\mathrm{z}; $\beta\equiv\frac{\mu_{J}}{hJ}\approx$0.6717\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}; $c\approx$113.8\text{,}\mathrm{k}\mathrm{H}\mathrm{z}; $d\approx$57.68\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}; $I=1$ is the total nuclear spin angular momentum in units of $\hbar$; $J=1$ is the total rotational angular momentum in units of $\hbar$; $\mu_{I}$ is the nuclear magnetic moment of a single nucleus; and $\mu_{J}$ is the magnetic moment due to molecular rotation. The first two terms describe the interaction of the nuclear and rotational magnetic moments with the external magnetic field, the third term describes the nuclear spin-rotational magnetic interaction Ramsey (1952); Kellogg et al. (1939, 1940), and the terms proportional to $d$ describe the magnetic spin-spin interaction of the two nuclei Ramsey (1952); Kellogg et al. (1939, 1940).

VI.2 Experiment and Observables

While there are many possible experimental protocols, we focus on the full interferometer mode used by Godsi et al. Godsi et al. (2017). The experiment is performed by initiating a continuous flux of oH2 molecules through the apparatus and measuring the current of the ionization detector while varying the first and second control fields ($B_{1}$ and $B_{2}$ in Figure 3).

In particular, $B_{1}$ is set to a specific value while $B_{2}$ is varied around the point ${-}B_{1}$ (i.e. about the spin-echo condition). In principle, $B_{2}$ could also be set to vary around ${+}B_{1}$, where spin echoes have also been observed Litvin et al. (2019), but we choose to vary $B_{2}$ about $-B_{1}$ to match the relevant experiment by Godsi et al. Godsi et al. (2017). This variation of the magnetic fields results in oscillatory curves of the detector current versus $B_{2}$, as shown in Figure 4 (a–c). These oscillations reflect the interference pattern that occurs when the various wavepackets recombine after passing through the final control field (see Section II). This interference pattern contains information about how the individual hyperfine states of the molecule interact with the sample surface.

The $x$ directed magnetic fields of a solenoid changes the energies of all of the $N_{\mathrm{R}}=9$ hyperfine states and induces all ${N_{\mathrm{R}}\choose 2}=36$ possible transitions. The frequencies of these transitions depend on the magnitude of the magnetic fields. By changing the magnitude of the second magnetic field, we are able to probe the rates of change of these transition frequencies with the magnetic field: the (generalized) gyromagnetic ratios $\gamma_{ij}(B)=\left|\frac{\mathrm{d}f_{ij}(B)}{\mathrm{d}B}\right|$, where $f_{ij}\equiv\frac{1}{h}\mathrm{\Delta}E_{ij}=\frac{E_{i}-E_{j}}{h}$, and $E_{i}$ is the energy of Ramsey state $i$ Godsi et al. (2017). The Fourier transforms of the oscillatory curves that give these gyromagnetic ratios are shown in Figure 4 (d–f). To obtain these results, we assumed that the surface normal of the sample lies in the scattering plane defined by the two branches and bisects the angle defined by the same two branches, such that $\alpha^{\prime}=\alpha=\gamma^{\prime}=\gamma=0$, $\beta=3\pi/8$ and $\beta^{\prime}=-5\pi/8$, where $\beta+\beta^{\prime}=-\pi/4=-\theta$ (see Eqn. 39 and Figure 3). Given this geometry and the axis definitions (Figure 3), the propagation matrices are defined with $+k_{R}$ in the first branch and $-k_{R}$ in the second. We also assume that the scattering transfer matrix is the identity matrix and is independent of energy, i.e. we assume for the present calculation that the only impact of scattering is the change of propagation direction, as modelled with rotation matrices (Section V.2).

The location of each feature in the spectra is reflective of a gyromagnetic ratio and is independent of the molecule-surface interactions, being only a function of the hyperfine energy level structure of oH2. The relative height of each feature, however, is dependent on the molecule-surface interactions, as exemplified in the experimental spectrum shown in Figure 4 (f). From Figure 4, one can see that integrating over the velocity distribution is important to produce the spin-echo effect and to bring the observed signal closer to experiment.

A different spectrum can be obtained for every possible value of $B_{1}$ and then combined to form a 2D map of the generalized gyromagnetic ratios and their contributing amplitudes as a function of $B_{1}$, as shown in Figure 5. This protocol is equivalent to observing the scattering of molecules with different internal hyperfine states as different values of the magnetic field in the first branch produce different superpositions of the hyperfine states. One can clearly see both the magnetic field-dependence of the gyromagnetic ratios, the impact of integrating over the velocity distribution, and the stark similarities and differences between the experimental and theory plots.

We now examine the sensitivity of the calculated signals to various changes in the scattering transfer matrices. Figures 6 and 7 demonstrate the impact of random variations of the scattering transfer matrix $\mathbf{\Sigma}$ on the oscillatory plots (for $B_{1}=$440\text{,}\mathrm{gauss}$$) and their spectra, respectively. For simplicity, we keep the matrix elements of $\mathbf{\Sigma}$ independent of energy.

The first row of each figure (labelled RP, Random Phases) reflects the impact of differing phases imparted to each hyperfine state after scattering. Specifically, $\mathbf{\Sigma}=\bigoplus_{i=1}^{9}e^{i\theta_{i}}$ is a diagonal unitary matrix whose nine phases $\theta_{i}$ are randomly chosen from a uniform distribution of width $2\pi$. Such a form of scattering would result from purely elastic scattering where the different hyperfine states probe the surface for different lengths of time (i.e. each state penetrates to a different depth or encounters a resonance with a different lifetime). Significant differences in the relative peak amplitudes can already be seen at this point, indicating that the calculated signal is sensitive to these phases.

The second row of each figure (labelled RDA, Random Diagonal Amplitudes) reflects the impact of differing state losses due to scattering. Specifically, $\mathbf{\Sigma}=\bigoplus_{i=1}^{9}A_{i}$ is a diagonal matrix whose diagonal elements are randomly chosen from a uniform distribution on the interval $\left[0,1\right)$. This form models the impact of different losses of each hyperfine state to different scattering directions, reactions with the surface, or adsorbtion to the surface. Again, significant changes are observed, indicating sensitivity to these features.

The third and fourth rows (respectively labelled ROM, Random Orthogonal Matrices, and RUM, Random Unitary Matrices) probe the impact of inelastic (projection $m_{I}m_{J}$-changing) scattering on the calculated signal. For the third row, $\mathbf{\Sigma}$ is an orthogonal matrix randomly drawn according to the Haar measure on $O(9)$, while $\mathbf{\Sigma}$ is a unitary matrix randomly drawn according to the Haar measure on $U(9)$ for the fourth row of each figure. Here, randomly drawing according to the Haar measure can be understood as analogous to drawing from the “uniform distribution” over the space of possible matrices Mezzadri (2007). The orthogonal matrices model inelastic scattering where no relative phase changes occur, while the unitary matrices model inelastic scattering where relative phase changes do occur. In both cases, there is no loss of total population during scattering. Clearly, the calculated signals are also sensitive to inelastic scattering events, both with and without relative phase changes. Finally, we can see that the number of peaks in the signal between 430 and 450 gauss varies as a function of the scattering matrix (compare the RDA and RUM plots in Figure 6, for example).

VII Comparison with a Semi-Classical Method

The present approach is fully quantum mechanical, while, in their Supplemental Material, Godsi et al. Godsi et al. (2017) have described a semi-classical method for calculating $P_{detection}$ that they used to model the propagation of oH2 in their molecular hyperfine interferometer (see Ref. Litvin et al. (2019) for the case of spin 1/2 particles). This semi-classical method treats the internal degrees of freedom of the molecules quantum mechanically and the centre-of-mass motion classically. As a result, the momentum changes induced by the magnetic field are ignored and every internal state is described as propagating at the initial velocity $v_{0}$ of the molecule. The internal degrees of freedom are treated by applying the time evolution operator for the time period $t_{i}\equiv\frac{L_{i}}{v_{0}}$ spent in each magnetic field of length $L_{i}$. That is, the propagation is calculated in the molecular reference frame with a time-dependent Hamiltonian. Here, we compare the results of the semi-classical and fully-quantum approaches for oH2.

We compare the two methods under conditions close to the original application of the semi-classical method to flux-detection measurements Godsi et al. (2017). We work with a field profile as shown in Figure 3, but with the second arm assumed to be of zero length and $B_{2}=0$. The field $B_{1}$ is varied. For the sake of comparison, we also set the state selector and detector fields in the transfer matrix method to $100\,000\text{\,}\mathrm{gauss}$ so that the basis changes performed by the transfer matrix method out of and into these regions match well the Clebsch-Gordon transformation from $\ket{m_{I}m_{J}}$ to $\ket{Fm}$ and its inverse, as used by the semi-classical method. Note that the off-diagonal elements of the full discontinuity matrix $\mathbf{K}$ (30) are still only ${\sim}10^{-5}$ at $100\,000\text{\,}\mathrm{gauss}$, such that their neglect still does not invalidate our fully quantum formalism at these large field strengths. We also retain the rotation from the first branch to the second branch and set the scattering matrix to $\mathbb{1}_{9}$. To maximize sensitivity of the comparison, we use only a single velocity when calculating $P_{detection}$ in both methods. All other parameters, including the state selector and detector relative state probabilities $\eta_{m_{I}m_{J}}$ and $\kappa_{m_{I}m_{J}}$, are as per Appendix D. This allows for a test that includes all incoming and outgoing states and their relative phases at experimentally relevant conditions.

The signals $P_{detection}\left(B_{1}\right)$ are calculated from $B_{1}$ to $B_{1}+$10\text{,}\mathrm{gauss}$$ for various values of $B_{1}$. We include 1500 datapoints in these ten-gauss intervals. The calculated signals are compared between the two methods by calculating their relative absolute difference at identical conditions. This produces a relative absolute difference at each of the 1500 magnetic field points. We then calculate the maximum, mean, and median relative absolute difference over the ten-gauss interval. Figure 8 shows how these maximum, mean, and median values vary as a function of $B_{1}$. At low fields, there is no significant dependence on the magnetic field and the relative absolute difference is below the expected experimental error. This lack of dependence on $B_{1}$ is possibly due to some residual numerical error present in the implementation of one or both methods, which masks any underlying field-dependence. At approximately $460\text{\,}\mathrm{gauss}$, however, the error begins to increase with the magnetic field until it saturates at approximately $46\,000\text{\,}\mathrm{gauss}$ at a relative absolute difference of approximately one. This increase in error as a function of magnetic field points to a systematic difference between the two methods.

To illustrate the difference between the two approaches in more detail, we plot the Fourier transforms of the calculated signals at various magnetic field strengths in Figure 9. At field strengths below $1000\text{\,}\mathrm{gauss}$, little difference is observed. At higher field strengths, the feature locations agree, while the Fourier amplitudes differ. The feature locations are determined by the eigenvalues of the Hamiltonian, identical in both methods, while the amplitudes are a function of the relative phases and amplitudes of the wavefunction components. These amplitudes and phases are expected to differ between the two methods at sufficiently high fields because of the approximations made in the semi-classical method.

In particular, the semi-classical method accounts for most of the relative phase and amplitude changes induced by the controlling magnetic fields. It does this by time-evolving the internal state vector for times that correspond to the time $t_{i}$ spent in each magnetic field by a molecule moving at its unchanged initial velocity. However, the semi-classical method ignores the small changes in the molecular velocity caused by the magnetic fields. These changes to the velocity modify the time spent in each magnetic field for each individual component of the internal state vector. Thus, $t_{i}$ should depend on the internal state $\ket{R}$. It is not immediately clear how to include these state-dependent velocity changes into the semi-classical method, however.

At low fields, these velocity changes and the dependence of $t_{i}$ on $\ket{R}$ are negligible and the fully quantum calculations agree with the semi-classical results to at least 0.1% for fields below $1000\text{\,}\mathrm{gauss}$. However, this agreement can only be expected to occur for surfaces that do not change between the surface-impact events of the spatially-separated wavepacket components (discussed in Section II). The maximum temporal separation between these impact events, caused by the velocity changes, varies from a few to several hundreds of picoseconds. Many surfaces do change on this timescale, as has been measured in several helium-3 spin echo experiments Jardine et al. (2009a, b, a, b, a, b); Hedgeland et al. (2016); Godsi et al. (2015); Jardine et al. (2009a, a, b). In other words, the semi-classical method cannot be used to probe the dynamics of surfaces, while the method presented in this manuscript opens the possibility to account for the surface dynamics with molecular scattering experiments.

VIII Conclusion

In this paper, we have developed a theoretical framework for simulating a surface-sensitive molecular hyperfine interferometer. The approach treats the interferometer as an effective one-dimensional system, accounting for the real experimental geometry by rotating the quantization axis of the hyperfine states at the scattering point. The time evolution of the molecular states is described fully coherently and accounts for the mixing of the hyperfine states and momentum changes induced by the magnetic fields in the experiment. The present approach is fully quantum mechanical and includes a full description of the internal-state-dependent spatial superpositions imposed on the molecular wavepackets by the controlling magnetic fields. This opens the possibility for a description of molecular scattering experiments that aim to probe surface dynamics on the picosecond to hundreds of picosecond time scale. To build the framework, we have derived and implemented a transfer matrix formalism that accounts for the internal (hyperfine) degrees of freedom of molecules and that allows for efficient computation of the experimental signal.

In the present work, the molecule-surface interaction is accounted for by a scattering transfer matrix (a transformed version of the scattering matrix) that is suitable for the description of experiments where the surface changes either much more slowly or much more quickly than the molecule-surface or wavepacket-surface interaction times (i.e., the molecule-surface scattering event does not involve energy transfer between the surface and the molecule). The extension to arbitrary surface dynamics (currently under investigation) requires a time-dependent scattering transfer matrix that reflects the underlying time-dependence of the molecule-surface interaction potential. Such a formalism would naturally incorporate energy transfers between the surface and the molecule during the scattering event. We have demonstrated, using the specific case of oH2, how the different features of the time-independent scattering transfer matrix, such as the phases of the diagonal elements, impact the experimental signal. In addition, we have shown that the experimental signal is sensitive to off-diagonal scattering matrix elements describing collisions that change the projection quantum numbers of the molecular hyperfine states without energy transfer between the molecule and the surface.

The present approach also sets the stage for solving the inverse scattering problem in molecular hyperfine interferometry by means of machine learning approaches, such as Bayesian optimization Vargas-Hernández et al. (2019); Krems (2019). For example, one can use the results of the transfer-matrix computations presented here to train Gaussian process models of the predicted experimental signal Krems (2019). The difference between the experimental observations and the results of the transfer-matrix computations can then be minimized by varying the scattering matrix elements, as described in our previous work Vargas-Hernández et al. (2019). The results of Bayesian optimization will determine the properties of the scattering matrix elements compatible with a given experimental measurement. These scattering matrix properties can then be used to gain physical insight into molecule-surface interactions and surface properties. They can also be used to test approximations used in ab initio calculations.

The formalism presented here is general to all closed-shell molecules and is flexible to describe various experimental setups. It can be used to explore various experimental protocols and evaluate their effectiveness at determining various molecule-surface interactions and surface properties. Thus, this paper provides the theoretical framework necessary to interpret a wide range of molecular hyperfine interferometry experiments, which are poised to apply molecular beam techniques to provide new information about molecule-surface interactions, surface morphologies, and surface dynamics.

Acknowledgements

This work is supported by NSERC of Canada and the Horizon 2020 Research and Innovation Programme grant 772228. We acknowledge useful discussions with Helen Chadwick, Geert Jan-Kroes, Ilya Litvin, and Tsofar Maniv and thank Oded Godsi for allowing us to use the source code of his implementation of the semi-classical method.

Appendix A Schrödinger Equation for Eigenstate Coefficients

Using Eqn. 19, the time-independent Schroödinger equation is

[TABLE]

where $\hat{H}$ is the total Hamiltonian (2) in the current region, $\hat{K}\equiv\frac{{\hat{p}}^{2}}{2m}$, we use $\hat{H}^{\text{R}}$ as a shorthand for $\hat{H}^{\text{R}}(\vec{B}_{\mathrm{loc}})$, and the last line was multiplied by $\bra{x_{0}R_{0}}$.

The different terms can be evaluated as

[TABLE]

where $\ket{k}$ is a momentum state with wavenumber $k$, and $m$ is the mass of the molecule. The additional factor of $\left(2\pi\right)^{-1}$ in Eqn. (45) comes from $\braket{x}{k}\equiv\left(2\pi\right)^{-\frac{1}{2}}e^{ikx}$. After inserting these three equations into Eqn. (44) and evaluating most of the sums, we obtain

[TABLE]

Noting that $k^{2}e^{ikx_{0}}=-\frac{\partial^{2}}{\partial x_{0}^{2}}e^{ikx_{0}}$ and $\int\frac{\mathrm{d}k}{2\pi}\leavevmode\nobreak\ e^{ik\left(x_{0}-x\right)}=\delta(x_{0}-x)$, we obtain Eqn. (20) after the relabelling $x_{0}\to x$.

Appendix B Coefficient Relations Across a Discontinuity

Since $\Phi_{R}^{E\tilde{R}}(x)\in C^{1}(x)$ for a specific value of $R$ and given Eqn. (26), we get the defining equations for the continuity of the wavefunction as

[TABLE]

where $S_{R^{-}R^{+}}^{*}\equiv\braket{R^{+}}{R^{-}}$, $k_{R^{\pm}}\equiv\frac{\sqrt{2m\left(E-E_{\text{R}^{\pm}}\right)}}{\hbar}$, and $E_{\text{R}^{\pm}}\equiv\braket{R^{\pm}}{\hat{H}^{\text{R}}\boldsymbol{\left(\right.}\vec{B}(0^{\pm})\boldsymbol{\left.\right)}}{R^{\pm}}$. There are $N_{R}$ such equations, one for each value of $R^{+}$.

Correspondingly, the defining equations for the continuity of the first derivative of the coefficients are

[TABLE]

Solving Eqns. (49) and (50) for the coefficients $A_{R^{+}}$ and $B_{R^{+}}$, we obtain Eqns. (27) and (28).

Appendix C Computational parameters used for the application to ortho-hydrogen

We take the mean velocity $v_{0}=$1436.14\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1} and the velocity spread to be $4\text{\,}\%\leavevmode\nobreak\ \mathrm{F}\mathrm{W}\mathrm{H}\mathrm{M}$. When performing the integral of Eqn. ([18](#S4.E18)), we take a $k$-space grid spacing $\Delta k=$1\text{\times}{10}^{4}\text{\,}{\mathrm{cm}}^{-1} and integrate from $-7\sigma_{k}$ to $+7\sigma_{k}$, where $\sigma_{k}$ is the Gaussian width in momentum space as defined in Section III. For the magnetic field profile and the angles between the two branches of the apparatus, see Figure 3. The relative probabilities used for the state selector probabilities $P_{R_{0}}$ and the detector coefficients $c_{R_{D}}$ are given in Table 1.

Where applicable, the parameters above were chosen to match those in the supplementary information of Godsi et al. Godsi et al. (2017), apart for the relative probabilities in Table 1. The relative probabilities in Table 1 were obtained from improved semi-classical calculations of the molecular propagation through the magnetic lens Godsi et al. (2017); Krüger et al. (2018).

Appendix D Computational parameters used for the comparison with the semi-classical method

We take the mean velocity $v_{0}=$1436.14\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}$$. The relative probabilities used for the state selector probabilities $P_{R_{0}}$ and the detector coefficients $c_{R_{D}}$ are given in Table 2. Where applicable, the parameters were chosen to match those in the supplementary information of Godsi et al. Godsi et al. (2017).