Phase-Sensitive Vector Terahertz Electrometry from Precision Spectroscopy of Molecular Ions

Abstract

This article proposes a new method for sensing THz waves that can allow electric field measurements traceable to the International System of Units and to the fundamental physical constants by using the comparison between precision measurements with cold trapped HD⁺ ions and accurate predictions of molecular ion theory. The approach exploits the lightshifts induced on the two-photon rovibrational transition at 55.9 THz by a THz wave around 1.3 THz, which is off-resonantly coupled to the HD⁺ fundamental rotational transition. First, the direction and the magnitude of the static magnetic field applied to the ion trap is calibrated using Zeeman spectroscopy of HD⁺. Then, a set of lightshifts are converted into the amplitudes and the phases of the THz electric field components in an orthogonal laboratory frame by exploiting the sensitivity of the lightshifts to the intensity, the polarization and the detuning of the THz wave to the HD⁺ energy levels. The THz electric field measurement uncertainties are estimated for quantum projection noise-limited molecular ion frequency measurements with the current accuracy of molecular ion theory. The method has the potential to improve the sensitivity and accuracy of electric field metrology and may be extended to THz magnetic fields and to optical fields.

1. Introduction

The sensitive detection of electromagnetic fields has a broad range of fundamental and technological applications, spanning from the search for physics beyond the Standard Model and tests of fundamental symmetries to time-keeping, navigation, modern communications, geophysics, and medical imaging. A goal of the metrology community was to make the measurements traceable to the International System of Units (SI) and to the fundamental physical constants. The measurement standards based on atoms and molecules were used for the determination of a number of physical quantities [1,2,3]—for example, the time (s), the length (m), and the mass (kg)—and of many fundamental constants [4], including the Rydberg constant, the fine structure constant, and various particle masses.

Atom-based techniques have been extended to the electric fields [5,6], enabling SI traceable measurements with fine spatial resolution. Microwave power measurement has been performed using the atomic candle method [7]. The sensitivity of the Rydberg states [8] and the well-known interactions with the electromagnetic fields were exploited using the electromagnetically induced transparency (EIT) and the Autler–Townes (AT) splitting phenomena [5,6,9] to probe electric fields oscillating from the radiofrequency [10] to the sub-terahertz domain [11]. The SI traceability of the electric field magnitude was ensured by the frequency measurement of the AT splitting, which was assumed as a linear ac-Stark effect depending on the Planck constant and the known dipole moment of the transition between the atomic Rydberg states. Comparing to the previous radiofrequency calibrations performed with standard electromagnetic fields and antenna probes [12,13], which display sensitivities at the 1 mV/(Hz^1/2.cm) level, and fractional accuracies in the range of 5–20%, the method using the atomic Rydberg transitions has many advantages for moderate to high-field measurements [6,14]. An improvement in accuracy is allowed by measurements with a fractional uncertainty better than 1% [5,15,16]. However, this method has difficulty to probe weak electric fields where the AT splitting is unresolved. Using the absorption spectrum, this method allowed detecting an electric field with the lowest magnitude of 8 µV/cm and a sensitivity of 30 µV/(Hz^1/2.cm) [5]. The improvement of the sensitivity at the level of a few µV/(Hz^1/2.cm) was achieved using homodyne detection [17] and frequency modulation [18], which allowed detecting a weak electric field of 1.8 µV/cm with matched filtering of the spectra. That corresponds to the current limit of sensitivity given by the photon shot noise on the detector used to record the spectra, which is three orders of magnitude worse than the quantum projection noise limit of the atomic sensor [17,18]. The surpassing of the photon shot noise limitation is nontrivial [19]. The measurements with Rydberg atomic transitions were extended to vector microwave electrometry with 0.5° angular resolution [20] and to a Rydberg-atom-based mixer for microwave phase shift measurement with an uncertainty better than 0.5 rad [21], offering also the opportunity to detect weak microwave signals [22].

This work proposes a new method for phase-sensitive vector electrometry based on the comparison of high-precision measurements with molecular ions with very accurate molecular theory predictions. The hydrogen molecular ions (HMI) are three-body molecular systems with hundreds of long-lived rotation–vibration energy levels in the ground electronic state. Starting from the mid-1960s, HMI were investigated with different experimental techniques (see for example [23,24,25,26,27,28,29,30]). The rotation–vibration spectroscopy data were addressed in the field of the metrology of particle masses [31]. Trapping HMI in a radiofrequency trap and sympathetically cooling with laser-cooled Be⁺ ions allowed measuring rovibrational lines with uncertainties at the 10⁻⁹ level [27,28,30]. New methods for precision measurements with trapped molecular ions were introduced recently [32,33]. Specifically for HMI, the sub-Doppler spectroscopy techniques, namely two-photon spectroscopy [34,35] and trapped ion cluster transverse excitation spectroscopy [36], may improve significantly the resolution and the accuracy. Recently, the fundamental rotational transition of the HD⁺ ion was probed with a fractional full-width at half-maximum (FWHM) linewidth of 3 × 10⁻¹² and a fractional uncertainty of 1.3 × 10⁻¹¹ [37]. In addition, the two-photon spectroscopy was demonstrated on the (v, L) = (0, 3)→(9, 3) transition with a fractional uncertainty of 2.9 × 10⁻¹² [38]. The ab initio calculations of physical parameters of HD⁺ using a set of fundamental constants reached a fractional accuracy of 7.6 × 10⁻¹² for vibrational transitions ignoring spin-structure effects [39] and a fractional accuracy at the 10⁻⁴ level or better for dipole moments of E1 transitions [40].

Currently, there is a lack of methods for calibrated measurements of the electric field of THz waves. A THz electric field may be characterized from the measurements of the ac-Stark shift (lightshift) that it induces on a two-photon rovibrational transition of HD⁺. The HD⁺ ions benefit from relatively high dipole moments of the transitions between adjacent rotational levels. These transitions are not allowed in the case of the homonuclear molecular ions such as H₂⁺. The THz wave is coupled off-resonantly to an E1 rotational transition from/to the energy levels addressed in two-photon spectroscopy. The THz electric field is characterized using the comparison between the experimental lightshift and the ab initio value of the lightshift calculated as a function of the CODATA2014 fundamental constants. A set of scalar, phase-less accurate optical frequency measurements may be translated into the amplitudes and the phases of the Cartesian components of the THz electric field, which points to the potential of this method for the full characterization of the electromagnetic fields. This article proposes converting two-photon infrared laser spectroscopy measurements with an expected fractional uncertainty at the 10⁻¹² level into THz electric field measurements for which the SI traceability is ensured by estimating the uncertainties at each stage of the experimental measurements and of the theoretical calculations.

2. Material and Methods

2.1. Theoretical Model of the HD⁺ Energy Levels

The vibration and the rotation of the HD⁺ ions define the most important scales of the energy levels in the electronic ground state. The rovibrational energy, which is calculated ab initio by neglecting the spin-structure effects, is the sum of the nonrelativistic Schrödinger energy with a series expansion of corrections that takes into account the relativistic, radiative, and nuclear-size-related contributions. The calculations performed up to the order α⁷ and including partial corrections of the order α⁸ (α is the fine structure constant) have a fractional accuracy estimated at the 10⁻¹² level [39].

The interactions between the spins of the electron ${\overrightarrow{\mathrm{S}}}_{\mathrm{e}}$, proton ${\overrightarrow{\mathrm{I}}}_{\mathrm{p}}$, and deuteron ${\overrightarrow{\mathrm{I}}}_{\mathrm{d}}$ that compose the HD⁺ ion with the rotational angular momentum $\overrightarrow{\mathrm{L}}$ lead to a hyperfine structure of energy levels [41]. The hyperfine energy levels are calculated ab initio, assuming the angular momentum coupling scheme $\overrightarrow{\mathrm{F}}={\overrightarrow{\mathrm{S}}}_{\mathrm{e}}+{\overrightarrow{\mathrm{I}}}_{\mathrm{p}};\overrightarrow{\mathrm{S}}=\overrightarrow{\mathrm{F}}+{\overrightarrow{\mathrm{I}}}_{\mathrm{d}};\overrightarrow{\mathrm{J}}=\overrightarrow{\mathrm{L}}+\overrightarrow{\mathrm{S}}$, by using a Breit–Pauli spin Hamiltonian in the nonrelativistic limit at the order α². The eigenvalues are labeled with the quantum numbers |vLFSJJ_z> for the vibration v, the squared angular momenta ${\Vert \overrightarrow{\mathrm{L}}\Vert }^2,{\Vert \overrightarrow{\mathrm{F}}\Vert }^2,{\Vert \overrightarrow{\mathrm{S}}\Vert }^2,{\Vert \overrightarrow{\mathrm{J}}\Vert }^2$, and the z-axis projection of $\overrightarrow{\mathrm{J}}$. Upon application of a small external static magnetic field defining the quantization axis, the HD⁺ magnetic states are split. The magnetic energy levels are calculated ab initio using a nonrelativistic Hamiltonian describing the spin interaction with the external magnetic field [42]. The value of the Zeeman shift of an energy level may be derived approximately as a quadratic dependence in function of the magnitude of the magnetic field. The energy of a magnetic state reads:

$$\mathrm{E}\left({\mathrm{v},\mathrm{L},\mathrm{F},\mathrm{S},\mathrm{J},\mathrm{J}}_{\mathrm{z}}\right)={\mathrm{E}}_{\mathrm{rv}}\left(\mathrm{v},\mathrm{L}\right)+{\mathrm{E}}_{\mathrm{hf}}\left(\mathrm{v},\mathrm{L},\mathrm{F},\mathrm{S},\mathrm{J}\right)+{\mathsf{\Delta }\mathrm{E}}_{\mathrm{z}}\left({\mathrm{v},\mathrm{L},\mathrm{F},\mathrm{S},\mathrm{J},\mathrm{J}}_{\mathrm{z}};\mathrm{B}\right),$$

(1)

by adding the rovibrational, hyperfine, and Zeeman energy contributions, respectively.

2.2. Two-Photon Spectroscopy for Sensing Electromagnetic Fields

This work proposes investigating the electromagnetic fields by Doppler-free spectroscopy of the HD⁺ ions in a radiofrequency trap that allows narrow linewidths and small systematic shifts. The parameters characterizing the electromagnetic fields are assumed with no time dependence during the experimental measurements. The proposed experimental setup (Figure 1a) builds on previous approaches demonstrated for HD⁺ spectroscopy [36,38]. Typically, ≈10² HD⁺ ions are stored together with ≈10³ Be⁺ ions that are laser-cooled with a 313 nm continuous-wave laser. The electrostatic interactions embed both ion species in a Coulomb crystal and sympathetically cool the HD⁺ ion motional degrees of freedom at the 10 mK level. The number of the HD⁺ ions is monitored through the change of the fluorescence of the Be⁺ ions at 313 nm upon excitation of the secular motion of the HD⁺ ions in the trap. A small static magnetic field with no spatial gradient, generated with three coil pairs in the Helmholtz configuration driven independently with three current sources, is applied to the ion trap.

Figure 1b indicates the rotation–vibration HD⁺ energy levels addressed in the sensing scheme. Two-photon spectroscopy is performed on the rovibrational transition (v, L) = (0, 0)→(2, 0) of HD⁺ using a stationary wave from an infrared laser tuned around 55.909 THz. The laser source for two-photon spectroscopy should have an ultranarrow linewidth, be tunable, and be continuously monitored against a frequency standard. These requirements may be met, for example, using the linewidth transfer technique approach [43], with a stabilized frequency comb and frequency down-conversion from near-infrared to the 5 µm spectral region. A feasible detection method is to dissociate the population transferred to the (v, L) = (2, 0) level with a 175 nm laser, as it was discussed in [44], without probing additionally a rotational transition. The population of the HD⁺ ions is initially distributed among the rotational levels of the ground vibrational state with a dependence determined by the thermal equilibrium with the blackbody radiation at room temperature. The blackbody radiation continuously recycles the HD⁺ ions population among these levels by driving electric dipole transitions, which are quantified here with the appropriate Einstein rate coefficients for spontaneous emission as well as stimulated absorption and emission. A set of rate equations, describing all transitions driven by the lasers and the blackbody radiation and taking into account the natural lifetimes of the relevant HD⁺ energy levels, allows characterizing the time dependence of the HD⁺ ions population during the resonance-enhanced multiphoton dissociation (REMPD) detection scheme assuming the hyperfine-free approximation. The lineshape of the two-photon resonance, probed with a 5362 nm laser, displays a full-width at half-maximum linewidth of 20 Hz for a two-photon transition rate of 10 s⁻¹, a dissociation rate of 200 s⁻¹, and an REMPD time of 10 s [44]. A similar approach based on REMPD allowed recently demonstrating the two-photon infrared spectroscopy of the (v, L) = (0, 3)→(9, 3) transition of HD⁺ [38]. The REMPD detection scheme proposed here estimates that a single datapoint may be recorded in 30 s, and the two-photon rovibrational line of HD⁺, which is used as frequency reference, may be defined within ten minutes. These timescales are comparable with the timescales of the recording of the rotational lines [37]. It is assumed that the laser can be referenced to the HD⁺ line with an uncertainty at a fraction of its linewidth within a ten-minute timescale.

The transitions between the spin states with a maximum total angular momentum J and extreme values of the projection J_z = ±J (stretched states) of the levels (v, L) = (0, 0) and (v, L) = (2, 0) are weakly split by the magnetic field [42]. The lightshifts of these transitions are used as probes for THz electric fields coupled off-resonantly to Zeeman subcomponents of the rotational transition (v, L) = (0, 0)→(0, 1) at 1.315 THz or (v, L) = (2, 0)→(2, 1) at 1.197 THz, respectively. The electric quadrupole shifts of the sublevels of the L = 0 states, induced by couplings with other L = 0 sublevels driven by a THz electric field with a spatial gradient, vanish [45]. The electric quadruple shifts of the sublevels of the L = 0 states due to off-resonant couplings to sublevels of the L = 2 states are neglected. In addition, the same experimental approach based on two-photon rovibrational spectroscopy of HD⁺ may be exploited for characterization of the magnetic field in the ion trap. Precision Zeeman spectroscopy with a 5196 nm laser of a sensitive subcomponent of the (v, L) = (0, 0)→(2, 2) two-photon transition, using a detection based on the dissociation of the (v, L) = (2, 2) level, may be exploited for determination of the magnitude of the magnetic field.

The fractional uncertainty of the frequency measurement of a transition is estimated with the Allan variance for the quantum projection noise limit:

$${\mathsf{\sigma }}_{\mathrm{y}}\left(\mathsf{\tau }\right)=\frac{1}{\mathsf{\pi }\mathrm{Q}\sqrt{{\mathrm{N}}_{\mathrm{ion}}}}\times \sqrt{\frac{{\mathrm{T}}_{\mathrm{c}}}{\mathsf{\tau }}},$$

(2)

which is expressed with the quality factor of the two-photon transition Q = f_2ph/Δf_HWHM in terms of the half-linewidth Δf_HWHM. A single measurement with N_ion ions at the two-photon resonance is performed during the cycle time T_c. The measurements are averaged during the interrogation time τ. For a single ion spectroscopy experiment with T_c = τ, where the linewidth is the inverse of the radiative lifetime of the excited energy level [46], the frequency uncertainty for (v, L) = (0, 0)→(2, 0) line is estimated at 2.49 Hz and for (v, L) = (0, 0)→(2, 2) line at 2.57 Hz.

This proposal is based on precision spectroscopy in an ion trap, frequency control of the spectroscopy laser, and measurements against a frequency standard, which are approaches that were previously developed for atomic ion clocks in frequency metrology institutes. The setup is a tabletop experiment involving specific lasers and optical components. Particularly, ultrastable near-infrared lasers and a stabilized frequency comb, referenced to a frequency standard, may be down-converted to the 5 µm spectral region and then amplified for two-photon spectroscopy or exploited with an enhancement cavity. The laser radiation at 175 nm may be obtained by frequency quadrupling with nonlinear crystals from a 700 nm amplified laser system.

2.3. Coordinate Frames and External Fields

The coordinate frames used here are represented in Figure 1a. The external fields applied to the HD⁺ ions are characterized relative to a Cartesian Laboratory Coordinate Frame locally fixed to the Earth’s surface $LCF\left({\overrightarrow{\mathrm{e}}}_{\mathrm{x}},{\overrightarrow{\mathrm{e}}}_{\mathrm{y}},{\overrightarrow{\mathrm{e}}}_{\mathrm{z}}\right)$, and its axes are pointing to East–North–Up. The Z-axis is oriented along the Earth’s ellipsoid normal direction. The magnetic field in the experimental setup is controlled with three pairs of magnetic coils. Each pair of opposite coils is driven with the same current. The coil pairs, driven independently by three stabilized current sources, define the Coil Coordinate Frame $CCF\left({\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{x}},{\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{y}},{\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{z}}\right)$, which is not necessarily orthogonal. The orientations of the CCF axes relative to the LCF are defined with the Euler angles: $\left({\mathsf{\alpha }}_{\mathrm{z}},\mathsf{\pi }/2+{\mathsf{\alpha }}_{\mathrm{y}}\right)$ for ${\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{x}}$,$\left(\mathsf{\pi }/2,\mathsf{\pi }/2-{\mathsf{\alpha }}_{\mathrm{x}}\right)$ for ${\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{y}}$, and $\left(0,0\right)$ for ${\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{z}}$, respectively. The nonorthogonality of the CCF axes is accounted here relative to the LCF frame with three small angles ${\mathsf{\alpha }}_{\mathrm{x}}{,\mathsf{\alpha }}_{\mathrm{y}}{,\mathsf{\alpha }}_{\mathrm{z}}<<1$. The magnetic field vector generated by the coil pairs is expressed as:

$$\overrightarrow{\mathrm{B}}={\mathrm{k}}_{\mathrm{x}}{\mathrm{I}}_1{\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{x}}+{\mathrm{k}}_{\mathrm{y}}{\mathrm{I}}_2{\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{y}}+{\mathrm{k}}_{\mathrm{z}}{\mathrm{I}}_3{\overrightarrow{\mathrm{e}}}_{\mathrm{c},\mathrm{z}},$$

(3)

as a function of the current-to-field parameters k_x,y,z, and the coil currents I_1,2,3. The conversion between the components of a vector in the orthogonal and nonorthogonal frames can be performed using a linear transformation [47].

The THz electric field vector, expressed with the complex amplitude ${\mathrm{E}}_{\mathrm{THz}}$, complex polarization vector $\hat{\mathsf{\epsilon }}$, and angular frequency $\mathsf{\omega }$, is decomposed further in three orthogonal linearly polarized components, which read:

$$\overrightarrow{\mathrm{E}}\left(\mathrm{t}\right)=\frac{{\mathrm{E}}_{\mathrm{THz}}\hat{\mathrm{e}}}{2}\times {\mathrm{e}}^{-\mathrm{i}\mathsf{\omega }\mathrm{t}}+\mathrm{c}.\mathrm{c}.={\displaystyle \sum _{\mathrm{j}=\left\{\mathrm{x},\mathrm{y},\mathrm{z}\right\}}\frac{{\mathrm{E}}_{\mathrm{j}}{\overrightarrow{\mathrm{e}}}_{\mathrm{j}}}{2}\times {\mathrm{e}}^{-\mathrm{i}\left(\mathsf{\omega }\mathrm{t}+{\mathsf{\phi }}_{\mathrm{j}}\right)}}+\mathrm{c}.\mathrm{c}.,$$

(4)

as a function of real positive amplitudes ${\mathrm{E}}_{\mathrm{j}}$ and phases ${\mathsf{\phi }}_{\mathrm{j}}$. This dependence is represented by the polarization ellipse of the THz wave, in terms of three field amplitudes and two phases, in which the third phase is arbitrarily fixed to zero ${\mathsf{\phi }}_z=0$. To describe the interaction with the radiation, it is suitable to define the Cartesian Molecular Ion Coordinate Frame $MICF\left({\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{x}},{\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{y}},{\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{z}}\right)$ such that the direction of ${\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{z}}$ is along the direction of the magnetic field that defines the quantization axis. The THz electric field vector is expressed in the MICF using the standard components:

$$\overrightarrow{\mathrm{E}}\left(\mathrm{t}\right)={\displaystyle \sum _{\mathrm{q}=\left\{-1,0,1\right\}}{\left(-1\right)}^{\mathrm{q}}\frac{{\mathrm{E}}_{-\mathrm{q}}{\overrightarrow{\mathrm{e}}}_{\mathrm{q}}}{2}\times {\mathrm{e}}^{-\mathrm{i}\left(\mathsf{\omega }\mathrm{t}+{\mathsf{\phi }}_{-\mathrm{q}}\right)}}+\mathrm{c}.\mathrm{c}.,$$

(5)

with real positive amplitudes ${\mathrm{E}}_0{,\mathrm{E}}_{\pm 1}$ and phases ${\mathsf{\phi }}_0,{\mathsf{\phi }}_{\pm 1}$ having linear or circular polarizations defined with ${\overrightarrow{\mathrm{e}}}_0={\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{z}},{\overrightarrow{\mathrm{e}}}_{\pm 1}=\mp \left({\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{x}}\pm {\mathrm{i}\overrightarrow{\mathrm{e}}}_{\mathrm{mi},\mathrm{y}}\right)/\sqrt{2}$. The orientation of the MICF relative to the $LCF$ can be changed by varying the currents in the coils. The relative orientation of the two coordinate frames is defined with the Euler angles $\left(\mathsf{\alpha },\mathsf{\beta },\mathsf{\chi }=0\right)$ (the third angle being fixed here arbitrarily to zero). Each standard component of the THz electric field in the MICF, denoted with ${\mathrm{E}}_{\mathsf{\pi },{\mathsf{\sigma }}^{\pm }}^{\left(\mathsf{\alpha },\mathsf{\beta }\right)}$, couples off-resonantly to the π or σ^± Zeeman subcomponents of the (v, L) = (v, 0)→(v, 1) transitions and induces a lightshift in proportion with its squared amplitude. The standard components of the THz electric field in the MICF are related to the polarization ellipse parameters in the LCF:

$$\begin{array}{ll}{\left({\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\alpha },\mathsf{\beta }\right)}\right)}^2& ={\mathrm{E}}_{\mathrm{z}}^2\cdot {\cos }^2\left(\mathsf{\beta }\right)+\left({\mathrm{E}}_{\mathrm{x}}^2\cdot {\cos }^2\left(\mathsf{\alpha }\right)+{\mathrm{E}}_{\mathrm{y}}^2\cdot {\sin }^2\left(\mathsf{\alpha }\right)\right)\cdot {\sin }^2\left(\mathsf{\beta }\right)\\ & +{\mathrm{E}}_{\mathrm{z}}\cdot {\mathrm{E}}_{\mathrm{x}}\cdot \sin \left(2\mathsf{\beta }\right)\cdot \cos \left(\mathsf{\alpha }\right)\cdot \cos \left({\mathsf{\phi }}_{\mathrm{x}}\right)+{\mathrm{E}}_{\mathrm{z}}\cdot {\mathrm{E}}_{\mathrm{y}}\cdot \sin \left(2\mathsf{\beta }\right)\cdot \sin \left(\mathsf{\alpha }\right)\cdot \cos \left({\mathsf{\phi }}_{\mathrm{y}}\right)\\ & +{\mathrm{E}}_{\mathrm{y}}\cdot {\mathrm{E}}_{\mathrm{x}}\cdot \sin \left(2\mathsf{\alpha }\right)\cdot {\sin }^2\left(\mathsf{\beta }\right)\cdot \cos \left({\mathsf{\phi }}_{\mathrm{x}}-{\mathsf{\phi }}_{\mathrm{y}}\right)\\ {\left({\mathrm{E}}_{{\mathsf{\sigma }}^{\pm }}^{\left(\mathsf{\alpha },\mathsf{\beta }\right)}\right)}^2& =\frac{1}{2}\left({\mathrm{E}}_{\mathrm{x}}^2+{\mathrm{E}}_{\mathrm{y}}^2+{\mathrm{E}}_{\mathrm{z}}^2-{\left({\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\alpha },\mathsf{\beta }\right)}\right)}^2\right)\pm {\mathrm{E}}_{\mathrm{x}}\cdot {\mathrm{E}}_{\mathrm{y}}\cdot \cos \left(\mathsf{\beta }\right)\cdot \sin \left({\mathsf{\phi }}_{\mathrm{x}}-{\mathsf{\phi }}_{\mathrm{y}}\right)\\ & \mp {\mathrm{E}}_{\mathrm{x}}\cdot {\mathrm{E}}_{\mathrm{z}}\cdot \sin \left(\mathsf{\alpha }\right)\cdot \sin \left(\mathsf{\beta }\right)\cdot \sin \left({\mathsf{\phi }}_{\mathrm{x}}\right)\pm {\mathrm{E}}_{\mathrm{y}}\cdot {\mathrm{E}}_{\mathrm{z}}\cdot \cos \left(\mathsf{\alpha }\right)\cdot \sin \left(\mathsf{\beta }\right)\cdot \sin \left({\mathsf{\phi }}_{\mathrm{y}}\right)\end{array}$$

(6)

2.4. Lightshifts in HD⁺ Spectroscopy

The coupling between the HD⁺ energy levels and the THz electric field is expressed in the electric dipole approximation with the interaction Hamiltonian:

$$\mathrm{V}=-\overrightarrow{\mathrm{d}}\cdot \overrightarrow{\mathrm{E}}\left(\mathrm{t}\right)=-\frac{1}{2}{\mathrm{E}}_{\mathrm{THz}}\left(\overrightarrow{\mathrm{d}}\cdot \hat{\mathsf{\epsilon }}\right)\times {\mathrm{e}}^{-\mathrm{i}\mathsf{\omega }\mathrm{t}}+\mathrm{H}.\mathrm{c}.,$$

(7)

where $\overrightarrow{\mathrm{d}}$ is the electric dipole operator in the laboratory frame. When the THz electric field is far from the resonance with another energy level, the lightshift of an energy level, calculated with the second-order perturbation theory [48], reads:

$${\mathsf{\delta }\mathrm{E}}_{\mathrm{n}}=-\frac{1}{4}{\left|{\mathrm{E}}_{\mathrm{THz}}\right|}^2{\displaystyle \sum _{\mathrm{r}}\mathrm{Re}\left[\frac{\langle \mathrm{n}|{\left(\overrightarrow{\mathrm{d}}\cdot \hat{\mathsf{\epsilon }}\right)}^{+}|\mathrm{r}\rangle \langle \mathrm{r}|\left(\overrightarrow{\mathrm{d}}\cdot \hat{\mathsf{\epsilon }}\right)|\mathrm{n}\rangle }{{\mathrm{E}}_{\mathrm{r}}-{\mathrm{E}}_{\mathrm{n}}-\hslash \mathsf{\omega }-\mathrm{i}\hslash \frac{{\mathsf{\gamma }}_{\mathrm{r}}+{\mathsf{\gamma }}_{\mathrm{n}}}{2}}+\frac{\langle \mathrm{n}|\left(\overrightarrow{\mathrm{d}}\cdot \hat{\mathsf{\epsilon }}\right)|\mathrm{r}\rangle \langle \mathrm{r}|{\left(\overrightarrow{\mathrm{d}}\cdot \mathsf{\epsilon } ^\right)}^{+}|\mathrm{n}\rangle }{{\mathrm{E}}_{\mathrm{r}}-{\mathrm{E}}_{\mathrm{n}}+\hslash \mathsf{\omega }+\mathrm{i}\hslash \frac{{\mathsf{\gamma }}_{\mathrm{r}}+{\mathsf{\gamma }}_{\mathrm{n}}}{2}}\right]},$$

(8)

in terms of the matrix elements of the dipole operator, the unperturbed energy levels E_r, E_n, and their decay rates ${\mathsf{\gamma }}_{\mathrm{r}}{,\mathsf{\gamma }}_{\mathrm{n}}$. The matrix elements of the dipole operator have been calculated ab initio in the Born–Oppenheimer approximation using nonrelativistic three-body wavefunctions [40]. The lightshift can be related to the reduced matrix elements of the dipole operator between rovibrational states, using the tensor formalism applied to the cyclic components d_q of the dipole operator. Equation (8) can be expanded with the squared modules of the standard components of the THz electric field:

$$\begin{array}{l}{ \mathsf{\delta }\mathrm{E}}_{\mathrm{n}}=-\frac{1}{4}{\displaystyle \sum _{\mathrm{q}=\left\{-1,0,1\right\}}{\left(-1\right)}^{\mathrm{q}}{\left|{\mathrm{E}}_{-\mathrm{q}}\right|}^2}{\mathsf{\alpha }}_{\mathrm{n},\mathrm{q}}\left(\mathsf{\omega }\right)\\ { \mathsf{\alpha }}_{\mathrm{n},\mathrm{q}}\left(\mathsf{\omega }\right)=\mathrm{Re}\left[\langle \mathrm{n}|{\mathrm{d}}_{-\mathrm{q}}\cdot \frac{1}{{\underset{\_}{\mathrm{H}}}^{*}-{\underset{\_}{\mathrm{E}}}_{\mathrm{n}}-\hslash \mathsf{\omega }}\cdot {\mathrm{d}}_{\mathrm{q}}|\mathrm{n}\rangle +\langle \mathrm{n}|{\mathrm{d}}_{-\mathrm{q}}\cdot \frac{1}{\underset{\_}{\mathrm{H}}-{\underset{\_}{\mathrm{E}}}_{\mathrm{n}}^{*}+\hslash \mathsf{\omega }}\cdot {\mathrm{d}}_{\mathrm{q}}|\mathrm{n}\rangle \right]\\ \end{array}$$

(9)

The last line of Equation (9) introduces the standard dynamic polarizabilities of the HD⁺ energy levels ${\mathsf{\alpha }}_{\mathrm{n},\mathrm{q}}\left(\mathsf{\omega }\right)$; these are expressed in terms of the complex operator $\underset{\_}{\mathrm{H}}=\mathrm{H}+\mathrm{i}\hslash \mathsf{\Gamma }$, which is defined as $\langle \mathrm{n}|\mathrm{H}|{\mathrm{n}}^{\prime }\rangle ={\mathsf{\delta }}_{{\mathrm{nn}}^{\prime }}{\mathrm{E}}_{\mathrm{n}}$ and $\langle \mathrm{n}|\mathsf{\Gamma }|{\mathrm{n}}^{\prime }\rangle ={\mathsf{\delta }}_{{\mathrm{nn}}^{\prime }}{\mathsf{\gamma }}_{\mathrm{n}}/2$, having the eigenvalues ${\underset{\_}{\mathrm{E}}}_{\mathrm{n}}={\mathrm{E}}_{\mathrm{n}}+\mathrm{i}\hslash {\mathsf{\gamma }}_{\mathrm{n}}/2$.

The standard dynamic polarizabilities of the Zeeman sublevels of the (v, L) = (0, 0) and (v, L) = (2, 0) states, depending on

$${ \mathsf{\alpha }}_{{\mathrm{v},\mathrm{L},\mathrm{F},\mathrm{S},\mathrm{J},\mathrm{J}}_{\mathrm{z}};\mathsf{\pi }/{\sigma }^{\pm 1}}\left({\mathrm{q},\mathrm{B},\mathrm{f}}_{\mathrm{THz}}\right)=\mathrm{f}\left(\left\{{\mathrm{E}}_{\mathrm{rv}}^{\left({\mathrm{v}}^{\prime },{\mathrm{L}}^{\prime }\right)}{,\mathrm{E}}_{\mathrm{hf}}^{\left({\mathrm{v}}^{\prime },{\mathrm{L}}^{\prime },{\mathrm{F}}^{\prime },{\mathrm{S}}^{\prime },{\mathrm{J}}^{\prime }\right)}{,\mathrm{Z}}_{\mathrm{k}}^{\left({\mathrm{v}}^{\prime },{\mathrm{L}}^{\prime },{\mathrm{F}}^{\prime },{\mathrm{S}}^{\prime },{\mathrm{J}}^{\prime }\right)}{,\mathsf{\mu },\mathsf{\gamma }}_{\left({\mathrm{v}}^{\prime },{\mathrm{L}}^{\prime }\right)}{,{\mathrm{J}}^{\prime }}_{\mathrm{z}}\right\}{;\mathrm{q},\mathrm{B},\mathrm{f}}_{\mathrm{THz}}\right),$$

(10)

are calculated approximately with Equation (9), by summing over the allowed electric dipole couplings to the sublevels of the (v, L) = (0, 1) and (v, L) = (2, 1) states, respectively. The standard dynamic polarizability is a function of a set $\left\{U_i\right\}$ of five types of theoretical parameters—the rovibrational energies (2 parameters from [49]), the hyperfine energies (up to 11 parameters from [41]), the Zeeman energies calculated approximately with quadratic magnetic field dependences with the parameters Z_{k = 1,2,3} (up to 33 parameters from [42]) and J_z, the reduced dipole moment (1 parameter from [40]), the natural linewidths of the energy levels (2 parameters from [46])—and a set of experimental parameters: the magnitude of the magnetic field B, the polarization q, and the frequency of the THz wave f_THz.

In order to avoid divergences and to maintain the approximation of a far-detuned THz electric field, the contribution of the resonant coupling was neglected on a small frequency domain centered on each resonance (assumed at 10 Hz for the polarizabilities of the (v, L) = (0, 0) sublevels, and at 1.1 kHz for the polarizabilities of the (v, L) = (2, 0) sublevels, respectively).

The covariances between the differential standard dynamic polarizabilities are calculated on the basis of the analytical dependences from Equation (10) by using the error propagation law, the uncertainties of the theoretical and experimental parameters, and their correlation coefficients. The uncertainties of the theoretical rovibrational energies are estimated at a fraction $\mathrm{u}\left({\mathrm{E}}_{\mathrm{rv}}\right)={10}^{-12}\cdot {\mathrm{E}}_{\mathrm{rv}}$ of their predicted values [49]. The uncertainty of the theoretical hyperfine energies from [41] is assumed at $\mathrm{u}\left({\mathrm{E}}_{\mathrm{hf}}\right)=500 \mathrm{Hz}$. The uncertainty of the theoretical Zeeman parameters from [42] is assumed at $\mathrm{u}\left({\mathrm{Z}}_{1,2}\right)={50 \mathrm{MHz}/\mathrm{T}}^2, \mathrm{u}\left({\mathrm{Z}}_3\right)=5 \mathrm{kHz}/\mathrm{T}$. The uncertainty of the theoretical dipole moments from [40] is assumed at $\mathrm{u}\left({\mathsf{\mu }}_{\mathrm{vL},{\mathrm{v}}^{\prime }{\mathrm{L}}^{\prime }}\right)=1.3\times {10}^{-4} \mathrm{a}.\mathrm{u}.$. The uncertainties of the theoretical radiative linewidths of the energy levels are estimated at $\mathrm{u}\left({\mathsf{\gamma }}_{\left(\mathrm{v},\mathrm{L}\right)=\left(0,1\right)}\right)=2.5\times {10}^{-7} \mathrm{Hz}$ and $\mathrm{u}\left({\mathsf{\gamma }}_{\left(\mathrm{v},\mathrm{L}\right)=\left(2,0\right)}\right)=\mathrm{u}\left({\mathsf{\gamma }}_{\left(\mathrm{v},\mathrm{L}\right)=\left(2,1\right)}\right)=0.49 \mathrm{Hz}$, by using values of energy level lifetimes from [46]. In addition, the correlation coefficients $\mathrm{corr}\left({\mathrm{E}}_{\mathrm{rv}}{,{\mathrm{E}}^{\prime }}_{\mathrm{rv}}\right)=\mathrm{corr}\left({\mathrm{E}}_{\mathrm{hf}}{,{\mathrm{E}}^{\prime }}_{\mathrm{hf}}\right)=\mathrm{corr}\left({\mathrm{Z}}_{{\mathrm{k}}_1}{,{\mathrm{Z}}^{\prime }}_{{\mathrm{k}}_2}\right)=\mathrm{corr}\left(\mathsf{\mu },{\mathsf{\mu }}^{\prime }\right)=\mathrm{corr}\left(\mathsf{\gamma },{\mathsf{\gamma }}^{\prime }\right)=1$ between these parameters are assumed to be equal to 1. The uncertainty of the magnitude of the magnetic field is calculated with the error propagation law by exploiting its dependence on the theoretical parameters and on the Zeeman-shifted molecular ion frequency (the approach is presented in Section 3.1). Finally, the uncertainty of the THz wave frequency is assumed at a fraction $\mathrm{u}\left({\mathrm{f}}_{\mathrm{THz}}\right)={10}^{-12}\cdot {\mathrm{f}}_{\mathrm{THz}}$ of the experimental value.

The standard dynamic polarizability of a stretched state of the ground rotational level and its uncertainty are calculated here for a given magnetic field as a function of the THz wave frequency and plotted in Figure 2. The profile of the standard dynamic polarizability has a narrow resonance when the THz wave is resonant with the transition (v, L, F, S, J, J_z) = (0, 0, 1, 2, 2, 2)→(0, 1, 1, 2, 3, 3). The uncertainty increases dramatically at resonance. Using a small detuning to the resonance for sensing a THz wave enhances the sensitivity at the expense of a loss in the calibration accuracy.

3. Results

3.1. Determination of the Magnetic Field Vector

The magnitude of the magnetic field in the ion trap can be determined from the measurements of the Zeeman shift of a two-photon rovibrational transition (v, L, F, S, J, J_z) = (n, J_z)→(n’, J’_z) of HD⁺ by adopting the linear approximation:

$$\begin{array}{l}\Vert \stackrel{\rightharpoonup }{\mathrm{B}}\left(\left\{{\mathrm{V}}_{\mathrm{k}}\right\}{;\mathrm{I}}_1{,\mathrm{I}}_2{,\mathrm{I}}_3\right)\Vert =\frac{{\mathrm{f}}_{\exp }-{\mathrm{f}}_{\mathrm{th}}}{{\mathsf{\eta }}_{\mathrm{th}}}\\ {\mathsf{\eta }}_{\mathrm{th}}=\left({\mathrm{Z}}_3^{\left({\mathrm{n}}^{\prime }\right)}{{\mathrm{J}}^{\prime }}_{\mathrm{z}}-{\mathrm{Z}}_3^{\left(\mathrm{n}\right)}{\mathrm{J}}_{\mathrm{z}}\right)/2\end{array}$$

(11)

as a function of the experimental frequency ${\mathrm{f}}_{\exp }$ that is shifted by the magnetic field from the theoretical ${\mathrm{f}}_{\mathrm{th}}$ zero-field hyperfine frequency. The linear Zeeman coefficient of the two-photon transition ${\mathsf{\eta }}_{\mathrm{th}}$ is expressed on the second line as a function of the theoretical parameters ${\mathrm{Z}}_3^{\left(\mathrm{n},{\mathrm{n}}^{\prime }\right)}$ for the Zeeman effect [42].

The magnetic field vector is controlled with three currents $\left({\mathrm{I}}_1{,\mathrm{I}}_2{,\mathrm{I}}_3\right)$ that drive the coil pairs independently. The uncertainty of the current in a coil pair is assumed at a fraction $\mathrm{u}\left({\mathrm{I}}_{1,2,3}\right)={10}^{-3}\left|{\mathrm{I}}_{1,2,3}\right|$ of the setting value. The magnetic field vector is fully described with a set $\left\{{\mathrm{V}}_{\mathrm{k}}\right\}$ of nine experimental parameters: the current-to-field parameters $\left({\mathrm{k}}_1{,\mathrm{k}}_2{,\mathrm{k}}_3\right),$ the nonorthogonality angles $\left({\mathsf{\alpha }}_{\mathrm{x}},{\mathsf{\alpha }}_{\mathrm{y}},{\mathsf{\alpha }}_{\mathrm{z}}\right)$ of the CCF, and the external bias magnetic field components $\left({\mathrm{B}}_{01}{,\mathrm{B}}_{02}{,\mathrm{B}}_{03}\right)$. As a function of these parameters, the magnetic field components in the LCF read:

$$\begin{array}{l}\stackrel{\rightharpoonup }{\mathrm{B}}\left(\left\{{\mathrm{k}}_1{,\mathrm{k}}_2{,\mathrm{k}}_3{,\mathsf{\alpha }}_{\mathrm{x}}{,\mathsf{\alpha }}_{\mathrm{y}}{,\mathsf{\alpha }}_{\mathrm{z}}{,\mathrm{B}}_{01}{,\mathrm{B}}_{02}{,\mathrm{B}}_{03}\right\}{;\mathrm{I}}_1{,\mathrm{I}}_2{,\mathrm{I}}_3\right)=\left({\mathrm{k}}_1{\mathrm{I}}_1-{\mathrm{k}}_2{\mathrm{I}}_2{\mathsf{\alpha }}_{\mathrm{z}}+{\mathrm{k}}_3{\mathrm{I}}_3{\mathsf{\alpha }}_{\mathrm{y}}+{\mathrm{B}}_{01}\right){\overrightarrow{\mathrm{e}}}_{\mathrm{x}}\\ +\left({\mathrm{k}}_2{\mathrm{I}}_2-{\mathrm{k}}_3{\mathrm{I}}_3{\mathsf{\alpha }}_{\mathrm{x}}+{\mathrm{B}}_{02}\right){\overrightarrow{\mathrm{e}}}_{\mathrm{y}}+\left({\mathrm{k}}_3{\mathrm{I}}_3+{\mathrm{B}}_{03}\right){\overrightarrow{\mathrm{e}}}_{\mathrm{z}}\end{array}$$

(12)

The set of experimental parameters can be calibrated by inversing Equation (11) using a nonlinear least-squares minimization. A number of N > 9 Zeeman spectroscopy measurements provide the Zeeman-shifted frequencies ${\mathrm{f}}_{\exp ,\mathrm{i}}$ at different current setpoints $\left({\mathrm{I}}_{1,\mathrm{i}}{,\mathrm{I}}_{2,\mathrm{i}}{,\mathrm{I}}_{3,\mathrm{i}}\right)$. The uncertainty of the calibration procedure is evaluated with the 9 × 9 covariance matrix ${\mathrm{G}}_{\mathrm{B}}$ of the estimation errors of the model parameters:

$${\mathrm{G}}_{\mathrm{B}}={\left({\mathrm{J}}_{\mathrm{B}}^{\mathrm{T}}\cdot {\mathrm{Y}}_{\mathrm{B}}^{-1}\cdot {\mathrm{J}}_{\mathrm{B}}\right)}^{-1},$$

(13)

which is expressed as a function of the N × N covariance matrix ${\mathrm{Y}}_{\mathrm{B}}$ of the input data ${\mathrm{y}}_{\mathrm{i}}=\left({\mathrm{f}}_{\exp ,\mathrm{i}}-{\mathrm{f}}_{\mathrm{th},\mathrm{i}}\right)/{\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}$ and the N × 9 Jacobian matrix ${\mathrm{J}}_{\mathrm{B}}=\left\{{{\left(\partial \Vert {\stackrel{\rightharpoonup }{\mathrm{B}}}_{\mathrm{i}}\Vert /\partial {\mathrm{V}}_{\mathrm{k}}\right)}_{\mathrm{i},\mathrm{k}}|}_{\left\{{\mathrm{V}}_{\mathrm{k}}\right\}=\left\{{\mathrm{V}}_{\mathrm{k};0}\right\}}\right\}$, which is calculated for the input data with the assumed values $\left\{{\mathrm{V}}_{\mathrm{k};0}\right\}$ of the model parameters. The value of the theoretical frequency at zero-magnetic field is expressed as the sum ${\mathrm{f}}_{\mathrm{th},\mathrm{i}}={\mathrm{f}}_{\mathrm{rv},\mathrm{i}}+{\mathrm{f}}_{\mathrm{hf},\mathrm{i}}$ of the rovibrational frequency and the hyperfine frequency. The errors of the input data are estimated as the quadratic sum of the contributions from the experimental uncertainty of the Zeeman subcomponent $\mathrm{u}\left({\mathrm{f}}_{\exp ,\mathrm{i}}\right)=2.57 \mathrm{Hz},$ and from the theoretical uncertainties of the rovibrational frequency $\mathrm{u}\left({\mathrm{f}}_{\mathrm{rv},\mathrm{i}}\right)={10}^{-12}\cdot {\mathrm{f}}_{\mathrm{rv},\mathrm{i}},$ the hyperfine frequency $\mathrm{u}\left({\mathrm{f}}_{\mathrm{hf},\mathrm{i}}\right)=0.5 \mathrm{kHz},$ and the Zeeman shift coefficient $\mathrm{u}\left({\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}\right)=5 \mathrm{kHz}/\mathrm{T}$. Nondiagonal elements in the covariance matrix of the input data arise from covariances between the theoretical parameters. The matrix ${\mathrm{Y}}_{\mathrm{B}}=\left\{{\mathrm{y}}_{\mathrm{b},\mathrm{ij}}\right\}$ reads:

$$\begin{array}{ll}{\mathrm{y}}_{\mathrm{b},\mathrm{ij}}& ={\left(\frac{\mathrm{u}\left({\mathrm{f}}_{\exp ,\mathrm{i}}\right)}{{\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}}\right)}^2\cdot {\mathsf{\delta }}_{\mathrm{ij}}+\frac{\mathrm{u}\left({\mathrm{f}}_{\mathrm{rv},\mathrm{i}}\right)\cdot \mathrm{u}\left({\mathrm{f}}_{\mathrm{rv},\mathrm{j}}\right)}{{\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}\cdot {\mathsf{\eta }}_{\mathrm{th},\mathrm{j}}}+\frac{\mathrm{u}\left({\mathrm{f}}_{\mathrm{hf},\mathrm{i}}\right)\cdot \mathrm{u}\left({\mathrm{f}}_{\mathrm{hf},\mathrm{j}}\right)}{{\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}\cdot {\mathsf{\eta }}_{\mathrm{th},\mathrm{j}}}\\ & +\left(\frac{\mathrm{u}\left({\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}\right)\cdot \left({\mathrm{f}}_{\mathrm{rv},\mathrm{i}}+{\mathrm{f}}_{\mathrm{hf},\mathrm{i}}\right)}{{\mathsf{\eta }}_{\mathrm{th},\mathrm{i}}^2}\right)\cdot \left(\frac{\mathrm{u}\left({\mathsf{\eta }}_{\mathrm{th},\mathrm{j}}\right)\cdot \left({\mathrm{f}}_{\mathrm{rv},\mathrm{j}}+{\mathrm{f}}_{\mathrm{hf},\mathrm{j}}\right)}{{\mathsf{\eta }}_{\mathrm{th},\mathrm{j}}^2}\right)\end{array}$$

(14)

Let us consider an experiment designed with the following current-to-field parameters ${\mathrm{k}}_1={\mathrm{k}}_2={\mathrm{k}}_3=-{10}^{-4}\mathrm{T}/\mathrm{A},$ which is operated under the Earth’s magnetic field with components along the directions east ${\mathrm{B}}_{01}=3.82\times {10}^{-7} \mathrm{T},$ north ${\mathrm{B}}_{02}=20.8733\times {10}^{-6} \mathrm{T},$ and up ${\mathrm{B}}_{03}=43.4494\times {10}^{-6} \mathrm{T},$ with the nonorthogonality angles of the CCF assumed as ${\mathsf{\alpha }}_{\mathrm{x}}={\mathsf{\alpha }}_{\mathrm{y}}={\mathsf{\alpha }}_{\mathrm{z}}=50 \mathrm{mrad}.$ The calibration is performed by Zeeman spectroscopy of the transition (v, L, F, S, J, J_z) = (0, 0, 1, 2, 2, −2)→(2, 2, 1, 2, 4, 0) of HD⁺. The frequencies of this Zeeman subcomponent are measured for 27 sets of current intensities $\left({\mathrm{I}}_1{,\mathrm{I}}_2{,\mathrm{I}}_3\right)$ = $\left({\mathrm{I}}_{\mathrm{off}1}+{\mathrm{n}}_1{\mathrm{I}}_0{,\mathrm{I}}_{\mathrm{off}2}+{\mathrm{n}}_2{\mathrm{I}}_0{,\mathrm{I}}_{\mathrm{off}3}+{\mathrm{n}}_3{\mathrm{I}}_0\right)$, where I₀ = 1 A, I_off1 = 3.82 mA, I_off2 = 208 mA, I_off3 = 422 mA, and n_1,2,3 = 0,1,2. These frequency-standard calibrated measurements may be performed during a timespan estimated by five hours. The errors of the calibrated parameters $\left\{{\mathrm{V}}_{\mathrm{k};\mathrm{cal}}\right\}$, determined from the nonlinear adjustment of the Zeeman spectroscopy data, are estimated by the covariance matrix ${\mathrm{G}}_{\mathrm{B}}$ calculated using Equations (11)–(14). Particularly, the estimations for the uncertainties are:

$$\begin{array}{l}\mathrm{u}\left({\mathrm{k}}_1\right)=1.9\times {10}^{-10} \mathrm{T}/\mathrm{A}, \mathrm{u}\left({\mathrm{k}}_2\right)=5.7\times {10}^{-10} \mathrm{T}/\mathrm{A}, \mathrm{u}\left({\mathrm{k}}_3\right)=1.1\times {10}^{-9} \mathrm{T}/\mathrm{A},\\ \mathrm{u}\left({\mathsf{\alpha }}_{\mathrm{x}}\right)=5.7\times {10}^{-6} \mathrm{rad}, \mathrm{u}\left({\mathsf{\alpha }}_{\mathrm{y}}\right)=4.0\times {10}^{-6} \mathrm{rad}, \mathrm{u}\left({\mathsf{\alpha }}_{\mathrm{z}}\right)=2.8\times {10}^{-6} \mathrm{rad},\\ \mathrm{u}\left({\mathrm{B}}_{01}\right)=4.7\times {10}^{-10} \mathrm{T}, \mathrm{u}\left({\mathrm{B}}_{02}\right)=4.4\times {10}^{-10} \mathrm{T}, \mathrm{u}\left({\mathrm{B}}_{03}\right)=5.6\times {10}^{-10} \mathrm{T},\end{array}$$

(15)

The Cartesian components $\mathrm{X}=\left({\mathrm{B}}_{\mathrm{x}}{,\mathrm{B}}_{\mathrm{y}}{,\mathrm{B}}_{\mathrm{z}}\right)$ of the magnetic field in the LCF or, alternatively, the components in spherical coordinates $\mathrm{X}=\left(\Vert \overrightarrow{\mathrm{B}}\Vert ,\mathsf{\alpha },\mathsf{\beta }\right)$ may subsequently be derived on the basis of Equation (12), using the calibrated parameters and the currents in the coil pairs. The angular components of the magnetic field define the Euler angles in the LCF for the direction of the quantization axis. The errors for these quantities are estimated with the covariances evaluated using the error propagation law:

$$\mathrm{cov}\left({\mathrm{X}}_{\mathrm{i}}\left(\left\{{\mathrm{V}}_{\mathrm{cal}}\right\}{,\mathrm{I}}_{1,\mathrm{i}}{,\mathrm{I}}_{2,\mathrm{i}}{\mathrm{I}}_{3,\mathrm{i}}\right){,\mathrm{X}}_{\mathrm{j}}\left(\left\{{\mathrm{V}}_{\mathrm{cal}}\right\}{,\mathrm{I}}_{1,\mathrm{j}}{,\mathrm{I}}_{2,\mathrm{j}}{\mathrm{I}}_{3,\mathrm{j}}\right)\right)={\displaystyle \sum _{\mathrm{p},\mathrm{q}=1}^9\frac{\partial {\mathrm{X}}_{\mathrm{i}}}{\partial {\mathrm{V}}_{\mathrm{p}}}\frac{\partial {\mathrm{X}}_{\mathrm{j}}}{\partial {\mathrm{V}}_{\mathrm{q}}}}{\mathrm{G}}_{\mathrm{B},\mathrm{pq}}+{\displaystyle \sum _{\mathrm{k}=1}^3\frac{\partial {\mathrm{X}}_{\mathrm{i}}}{\partial {\mathrm{I}}_{\mathrm{k}}}\frac{\partial {\mathrm{X}}_{\mathrm{j}}}{\partial {\mathrm{I}}_{\mathrm{k}}}}\mathrm{u}\left({\mathrm{I}}_{\mathrm{k},\mathrm{i}}\right)\mathrm{u}\left({\mathrm{I}}_{\mathrm{k},\mathrm{j}}\right),$$

(16)

The previous equation expresses the covariance as the sum of the contributions from the experimental parameters determined with the magnetic field calibration procedure and from setting currents in the coil pairs.

The Earth’s magnetic field may be cancelled in the ion trap by setting the current intensities that yield a null result in Equation (12): I_null1 = −6.34 mA, I_null2 = 230 mA, and I_null3 = 434 mA. The uncertainty is estimated at 85 nT, with the root sum of squares of the uncertainties of the Cartesian components of the magnetic field in the LCF. The contribution of the Earth’s magnetic field to the total magnetic field components in the LCF is canceled when the total current in each coil pair is expressed as the sum of the nulling current with an offset current. The total currents read $\left({\mathrm{I}}_{\mathrm{null}1}+{\mathrm{I}}_{\mathrm{offset}1}{,\mathrm{I}}_{\mathrm{null}2}+{\mathrm{I}}_{\mathrm{offset}2}{,\mathrm{I}}_{\mathrm{null}3}+{\mathrm{I}}_{\mathrm{offset}3}\right)$ = $\left({\mathrm{I}}_{\mathrm{null}1}+{\mathrm{I}}_0\cdot \sin \mathsf{\theta }\cdot \cos \mathsf{\varphi },{ \mathrm{I}}_{\mathrm{null}2}+{\mathrm{I}}_0\cdot \sin \mathsf{\theta }\cdot \sin \mathsf{\varphi },{ \mathrm{I}}_{\mathrm{null}3}+{\mathrm{I}}_0\cdot \cos \mathsf{\theta }\right)$ in Cartesian coordinate parametrization and in spherical coordinate parametrization, respectively.

The total uncertainties of the spherical components of the magnetic field in the LCF are calculated using Equations (12)–(16) and shown in Figure 3, in the case where the total currents in the coil pairs are described by the spherical coordinate parametrization with I₀ = 1 A. The uncertainties from the driving currents yield the dominant contribution to the uncertainty of the magnetic field components that is at the $\mathrm{u}\left(\Vert \overrightarrow{\mathrm{B}}\Vert ,\mathsf{\alpha },\mathsf{\beta }\right)=\left({10}^{-7} \mathrm{T}, {10}^{-3} \mathrm{rad}, {10}^{-3} \mathrm{rad}\right)$ level or better. The contribution to this uncertainty coming only from the experimental measurements, which is at the $\mathrm{u}\left(\Vert \overrightarrow{\mathrm{B}}\Vert ,\mathsf{\alpha },\mathsf{\beta }\right)=\left({10}^{-9} \mathrm{T}, {10}^{-6} \mathrm{rad}, {10}^{-5} \mathrm{rad}\right)$ level, is orders of magnitude smaller. Although the Euler angles do not depend on the value of the parameter I₀, their uncertainties do. As for the magnitude of the magnetic field, these uncertainties decrease by increasing the value of I₀.

3.2. Determination of the Polarization Ellipse of the THz Wave

The lightshifts induced on a two-photon rovibrational transition of HD⁺ by coupling the THz wave to the molecular ions ${\mathsf{\delta }\mathrm{f}}_{\mathrm{k}}^{\left(\mathsf{\alpha },\mathsf{\beta }\right)}={\mathrm{f}}_{\exp ,\mathrm{k}}^{\mathrm{THz} \mathrm{ON}}-{\mathrm{f}}_{\exp ,\mathrm{k}}^{\mathrm{THz} \mathrm{OFF}}$ are exploited here for characterization of the THz wave. For a given orientation of the magnetic field, three independent lightshift measurements may allow the determination of the squared modules of the standard components of the THz electric field in the MICF by solving a system of three equations depending on differential standard dynamic polarizabilities. The measurements should be performed with nondegenerate transitions that are experimentally resolved. In order to determine the amplitudes and phases $\left(E_x,E_y,E_z,{\mathsf{\phi }}_x,{\mathsf{\phi }}_y\right)$ of the THz electric field components in the LCF, a set of five independent measurements of the squared modules of the standard components of the electric field in the MICF $E_{\mathsf{\pi }}^{\left({\mathsf{\alpha }}_i,{\mathsf{\beta }}_i\right)},E_{{\mathsf{\sigma }}^{\pm }}^{\left({\mathsf{\alpha }}_i,{\mathsf{\beta }}_i\right)}$ are required for the inversion of Equation (6), using at least two different orientations of the magnetic field. The independence of the measurements means that the relevant systems of equations should be nonsingular. On one hand, that is possible by using different magnitudes of the magnetic field or by probing lightshifts on different two-photon hyperfine or Zeeman transitions; on the other hand, that is possible by using a set of suitable orientations of the magnetic field that takes into account the periodicity of the trigonometric functions in Equation (6). The contributions from offset THz electric fields in the LCF, similar to the blackbody radiation field, are not addressed here, and it is assumed that these fields do not vary when the THz wave is coupled to the molecular ions.

The lightshift measurements are performed by adjusting the magnitude of the magnetic field so that the THz wave is near resonant with a Zeeman subcomponent of the (v, L) = (v, 0)→(v, 1) dipole-allowed transitions of HD⁺. These resonance frequencies may be tuned up to ±7 GHz by applying a magnetic field up to 10⁻⁴ T. Using a small detuning to the resonance frequency is interesting for increasing the lightshift. The drawbacks are significant mixings of the Zeeman energy levels that may invalidate the second-order perturbation theory approach, a change of their populations by driving rotational transitions with the THz electric field, and a change of the quantization axis relative to the magnetic field orientation.

This work proposes exploiting six independent lightshift measurements ${\mathrm{F}}_{\mathrm{LS}}={\left\{{\mathsf{\delta }\mathrm{f}}_{\mathrm{i}}^{\left({\mathsf{\alpha }}_{\mathrm{i}}{,\mathsf{\beta }}_{\mathrm{i}}\right)}\right\}}_{\mathrm{i}=1,2,3,4,5,6}$ performed on a set of two-photon transitions:

$${ \mathsf{\delta }\mathrm{f}}_{\mathrm{i}}\left({\mathrm{I}}_{1,\mathrm{i}}{,\mathrm{I}}_{2,\mathrm{i}}{,\mathrm{I}}_{3,\mathrm{i}}\right)=-\frac{1}{4}{\displaystyle \sum _{\mathrm{q}=\left\{-1,0,1\right\}}{\left(-1\right)}^{\mathrm{q}}{\left|{\mathrm{E}}_{\mathrm{q}}^{\left({\mathsf{\alpha }}_{\mathrm{i}}{,\mathsf{\beta }}_{\mathrm{i}}\right)}\right|}^2}{\mathsf{\Delta }\mathsf{\alpha }}_{\mathrm{i},-\mathrm{q}}\left(\left\{{\mathrm{U}}_{\mathrm{j}}\right\},\Vert \overrightarrow{\mathrm{B}}\left(\left\{{\mathrm{V}}_{\mathrm{k}}\right\}{,\mathrm{I}}_{1,\mathrm{i}}{,\mathrm{I}}_{2,\mathrm{i}}{,\mathrm{I}}_{3,\mathrm{i}}\right)\Vert {,\mathrm{f}}_{\mathrm{THz}}\right),$$

(17)

by setting suitable currents $\left({\mathrm{I}}_{1,\mathrm{i}}{,\mathrm{I}}_{2,\mathrm{i}}{,\mathrm{I}}_{3,\mathrm{i}}\right)$ in the coil pairs such as $\left({\mathsf{\alpha }}_i,{\mathsf{\beta }}_i\right)=\left(0,\mathsf{\pi }/2\right)$ for i = 1, 2, 3 and $\left({\mathsf{\alpha }}_i,{\mathsf{\beta }}_i\right)=\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)$ for i = 4, 5, 6. These measurements may be performed during a timespan by two hours. These lightshifts are related to the squared modules of the standard components of the THz electric field, which are expressed as ${\mathrm{E}}_{\mathrm{STD}}^2=\left\{{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2, {\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2, {\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2, {\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2, {\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2, {\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2\right\},$ by using the corresponding differential standard dynamic polarizabilities ${\mathsf{\Delta }\mathsf{\alpha }}_{\mathrm{i},-\mathrm{q}}$ for the two-photon transitions. These are calculated as a function of the relevant set of theoretical parameters $\left\{{\mathrm{U}}_{\mathrm{j}}\right\}$, the magnitude of the magnetic field expressed with the experimental parameters {V_k} that was calibrated previously by Zeeman spectroscopy, and the frequency of the THz wave. Equation (17) can be written in the matrix form by introducing a 6 × 6 matrix $\mathrm{A}\left({\mathsf{\Delta }\mathsf{\alpha }}_{\mathrm{i},-\mathrm{q}}\right)$, which is composed of two 3 × 3 diagonal blocks having as matrix elements the differential standard dynamic polarizabilities multiplied with the suitable prefactors. The standard components of the THz electric field are derived by inversion ${\mathrm{E}}_{\mathrm{STD}}^2={\left[\mathrm{A}\left({\mathsf{\Delta }\mathsf{\alpha }}_{\mathrm{i},-\mathrm{q}}\right)\right]}^{-1}\cdot {\mathrm{F}}_{\mathrm{LS}}$, and their covariances are estimated from the covariances of the inverse matrix elements and the covariances of ${\mathrm{F}}_{\mathrm{LS}}$:

$$\begin{array}{l} \mathrm{cov}\left({\mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2{,\mathrm{E}}_{\mathrm{STD},\mathrm{j}}^2\right)={\displaystyle \sum _{\mathrm{k},\mathrm{l}}{\mathrm{A}}_{\mathrm{ik}}^{-1}{\mathrm{A}}_{\mathrm{jl}}^{-1}\mathrm{cov}\left({\mathrm{F}}_{\mathrm{LS},\mathrm{k}}{,\mathrm{F}}_{\mathrm{LS},\mathrm{l}}\right)}+{\displaystyle \sum _{\mathrm{s},\mathrm{t}}{\mathrm{F}}_{\mathrm{LS},\mathrm{s}}{\mathrm{F}}_{\mathrm{LS},\mathrm{t}}\mathrm{cov}\left({\mathrm{A}}_{\mathrm{is}}^{-1}{,\mathrm{A}}_{\mathrm{jt}}^{-1}\right)}\\ \mathrm{cov}\left({\mathrm{A}}_{{\mathrm{i}}_1{\mathrm{j}}_1}^{-1}{,\mathrm{A}}_{{\mathrm{i}}_2{\mathrm{j}}_2}^{-1}\right)={\displaystyle \sum _{\mathrm{p},\mathrm{r},\mathrm{s},\mathrm{t}}{\mathrm{A}}_{{\mathrm{i}}_1\mathrm{p}}^{-1}{\mathrm{A}}_{{\mathrm{rj}}_1}^{-1}{\mathrm{A}}_{{\mathrm{i}}_2\mathrm{s}}^{-1}{\mathrm{A}}_{{\mathrm{tj}}_2}^{-1}}\mathrm{cov}\left({\mathrm{A}}_{\mathrm{pr}}{,\mathrm{A}}_{\mathrm{st}}\right)\\ \mathrm{cov}\left({\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_1{,\mathrm{q}}_1}\left(\left(\left\{{\mathrm{U}}_{{\mathrm{i}}_1}\right\},\Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_1}\Vert {,\mathrm{f}}_{\mathrm{THz}}\right)\right){,\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_2{,\mathrm{q}}_2}\left(\left\{{\mathrm{U}}_{{\mathrm{i}}_2}\right\},\Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_2}\Vert {,\mathrm{f}}_{\mathrm{THz}}\right)\right)={\displaystyle \sum _{{\mathrm{i}}_1{,\mathrm{i}}_2}\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_1{,\mathrm{q}}_1}}{\partial {\mathrm{U}}_{{\mathrm{i}}_1}}}\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_2{,\mathrm{q}}_2}}{\partial {\mathrm{U}}_{{\mathrm{i}}_2}}\mathrm{cov}\left({\mathrm{U}}_{{\mathrm{i}}_1}{,\mathrm{U}}_{{\mathrm{i}}_2}\right)\\ +\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_1{,\mathrm{q}}_1}}{\partial \Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_1}\Vert }\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_2{,\mathrm{q}}_2}}{\partial \Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_2}\Vert }\mathrm{cov}\left(\Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_1}\Vert ,\Vert {\overrightarrow{\mathrm{B}}}_{{\mathrm{p}}_2}\Vert \right)+\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_1{,\mathrm{q}}_1}}{\partial {\mathrm{f}}_{\mathrm{THz}}}\frac{\partial {\mathsf{\Delta }\mathsf{\alpha }}_{{\mathrm{p}}_2{,\mathrm{q}}_2}}{\partial {\mathrm{f}}_{\mathrm{THz}}}{\mathrm{u}}^2\left({\mathrm{f}}_{\mathrm{THz}}\right)\end{array}$$

(18)

The errors of ${\mathrm{F}}_{\mathrm{LS}}$ are assumed to be uncorrelated, and the uncertainties $\mathrm{u}\left({\mathsf{\delta }\mathrm{f}}_{\mathrm{i}}\right)=\sqrt{2}\cdot 2.49 \mathrm{Hz}$ are the same. The second line of Equation (18) indicates the expression of the covariance of the inverse matrix elements calculated with the approach presented in [50]. The last two lines give the expression of the covariance between the differential standard dynamic polarizabilities, which are determined with the error propagation law from the covariances between the theoretical parameters $\left\{{\mathrm{U}}_{\mathrm{j}}\right\}$, the dependences on the magnetic field magnitude, and the frequency of the THz wave, respectively. The covariances of $\left\{{\mathrm{U}}_{\mathrm{j}}\right\}$ are indicated in Section 2.4, the covariance of the magnetic field magnitudes may be calculated with Equation (16), and the uncertainty of the THz wave frequency $\mathrm{u}\left({\mathrm{f}}_{\mathrm{THz}}\right)={10}^{-12}\cdot {\mathrm{f}}_{\mathrm{THz}}$ is assumed fractionally at the ppt level.

Furthermore, the amplitudes and phases of the THz electric field components in the LCF are derived from these standard components by inverting the nonlinear system of equations of Equation (6). The solutions ${\mathrm{X}}_0=\left({\mathrm{E}}_{\mathrm{x},0}{,\mathrm{E}}_{\mathrm{y},0}{,\mathrm{E}}_{\mathrm{z},0},{\mathsf{\phi }}_{\mathrm{x},0},{\mathsf{\phi }}_{\mathrm{y},0}\right)$ are expressed with the analytic dependences on ${\mathrm{E}}_{\mathrm{STD},0}^2=\left({\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2,{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2,{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2,{\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2,{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2\right)$ for this particular choice of orientations:

$$\begin{array}{l}{\mathrm{E}}_{\mathrm{x},0}=\sqrt{\left({\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2+{\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2+{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2\right)-\left({\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2+{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2\right)}\\ {\mathrm{E}}_{\mathrm{y},0}=\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|\\ {\mathrm{E}}_{\mathrm{z},0}=\sqrt{\left({\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2+{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2\right)-{\left|{\mathrm{E}}_{\mathsf{\pi }}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2}\\ {\mathsf{\phi }}_{\mathrm{x},0}=\arcsin \left[\frac{{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2-{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)}\right|}^2}{{2\mathrm{E}}_{\mathrm{x},0}{\mathrm{E}}_{\mathrm{z},0}}\right]\\ {\mathsf{\phi }}_{\mathrm{y},0}=\arcsin \left[\frac{{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{+}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2-{\left|{\mathrm{E}}_{{\mathsf{\sigma }}^{-}}^{\left(0,\mathsf{\pi }/2\right)}\right|}^2}{{2\mathrm{E}}_{\mathrm{y},0}{\mathrm{E}}_{\mathrm{z},0}}\right]\end{array}$$

(19)

In order to estimate the uncertainties, Equation (6) is linearized around these solutions:

$$\begin{array}{l}{\mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2={\mathrm{E}}_{\mathrm{STD},\mathrm{i},0}^2+{\displaystyle \sum _{\mathrm{k}=\left\{\mathrm{x},\mathrm{y},\mathrm{z}\right\}}\frac{\partial {\mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2}{\partial {\mathrm{E}}_{\mathrm{k}}}\left({\mathrm{E}}_{\mathrm{k}}-{\mathrm{E}}_{\mathrm{k},0}\right)}+{\displaystyle \sum _{\mathrm{k}=\left\{\mathrm{x},\mathrm{y}\right\}}\frac{\partial {\mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2}{\partial {\mathsf{\phi }}_{\mathrm{k}}}\left({\mathsf{\phi }}_{\mathrm{k}}-{\mathsf{\phi }}_{\mathrm{k},0}\right)}\\ { \mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2-{\mathrm{E}}_{\mathrm{STD},\mathrm{i},0}^2={\displaystyle \sum _{\mathrm{k}=1}^5{\left({\mathrm{J}}_{{\mathrm{E}}_{\mathrm{STD}}^2}\right)}_{\mathrm{ik}}\left({\mathrm{X}}_{\mathrm{k}}-{\mathrm{X}}_{\mathrm{k},0}\right)}\end{array}$$

(20)

The second line expresses the linear dependence in a compact form with the 6 × 6 Jacobian matrix ${\mathrm{J}}_{{\mathrm{E}}_{\mathrm{STD}}^2}=\left\{{{\left(\partial {\mathrm{E}}_{\mathrm{STD},\mathrm{i}}^2/\partial {\mathrm{X}}_{\mathrm{k}}\right)}_{\mathrm{i},\mathrm{k}}|}_{\left\{{\mathrm{X}}_{\mathrm{k}}\right\}=\left\{{\mathrm{X}}_{\mathrm{k},0}\right\}}\right\}$, which is calculated using the solution ${\mathrm{X}}_0$. The errors for the amplitudes and phases of the electric field components in the LCF are given by the 5 × 5 covariance matrix ${\mathrm{G}}_{\mathrm{E}}$ of the estimation errors of the model parameters:

$${\mathrm{G}}_{\mathrm{E}}={\left({\mathrm{J}}_{{\mathrm{E}}_{\mathrm{STD}}^2}^{\mathrm{T}}\cdot {\mathrm{Y}}_{{\mathrm{E}}_{\mathrm{STD}}^2}^{-1}\cdot {\mathrm{J}}_{{\mathrm{E}}_{\mathrm{STD}}^2}\right)}^{-1},$$

(21)

as a function of the 6 × 6 covariance matrix ${\mathrm{Y}}_{{\mathrm{E}}_{\mathrm{STD}}^2}$ of the standard components at $\left(\mathsf{\alpha },\mathsf{\beta }\right)=\left(0,\mathsf{\pi }/2\right)$ and $\left(\mathsf{\alpha },\mathsf{\beta }\right)=\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)$, which is calculated using Equation (18).

Let us discuss the application of this method for the characterization of a THz wave that is linearly polarized along an arbitrary direction. The amplitudes of the Cartesian components of the electric field are parametrized in the LCF as $\left({\mathrm{E}}_{\mathrm{x},0}{,\mathrm{E}}_{\mathrm{y},0}{,\mathrm{E}}_{\mathrm{z},0}\right)=\left({\mathrm{E}}_0\cdot \sin \mathsf{\theta }\cdot \cos \mathsf{\varphi },{ \mathrm{E}}_0\cdot \sin \mathsf{\theta }\cdot \sin \mathsf{\varphi },{ \mathrm{E}}_0\cdot \cos \mathsf{\theta }\right)$, as a function of the amplitude of the THz wave and two spherical angles $\left({\mathrm{E}}_0,\mathsf{\varphi },\mathsf{\theta }\right)$. The relevant phases are assumed as ${\mathsf{\phi }}_{\mathrm{x},0}={\mathsf{\phi }}_{\mathrm{y},0}=0$. Here, the frequency of the THz wave is assumed at f_THz = 1,314,947,502.3 kHz (a detuning of 1.6 MHz to the frequency of the (v, L, F, S, J) = (0, 0, 1, 2, 2)→(0, 1, 1, 2, 3) hyperfine transition at zero-magnetic field) and the intensity at 1 W/m² (E₀ = 27.42 mV/m). A set of six lightshifts are measured for the two-photon transition between the stretched states (v, L, F, S, J, J_z) = (0, 0, 1, 2, 2, 2)→(2, 0, 1, 2, 2, 2) of HD⁺ for two orientations of the magnetic field $\left(\mathsf{\alpha },\mathsf{\beta }\right)=\left(0,\mathsf{\pi }/2\right),\left(\mathsf{\alpha },\mathsf{\beta }\right)=\left(\mathsf{\pi }/2,\mathsf{\pi }/2\right)$ in the case of three different values of the magnitude of the magnetic field B₁ = 10⁻⁶ T, B₂ = 5 × 10⁻⁶ T, and B₃ = 10⁻⁵ T. The total uncertainties of the Cartesian components of the THz wave in the LCF and of their phases are calculated using Equations (16)–(21) on the basis of the standard dynamic polarizabilities given by Equation (10) and by using the calculated covariance matrix ${\mathrm{G}}_{\mathrm{B}}$ of the magnetic field components derived from Equations (11)–(14). The results are plotted in Figure 4 using the spherical angle parametrization $\left({\mathsf{\varphi }}_{\mathrm{i}},{\mathsf{\theta }}_{\mathrm{j}}\right)=\left({10}^{-11}+\left(\mathsf{\pi }\mathrm{i}/60\right){ \mathrm{rad}, 10}^{-11}+\left(\mathsf{\pi }\mathrm{j}/60\right) \mathrm{rad}\right)$ for i = 0–30, j = 0–30. The total uncertainties are estimated at:

$$\begin{array}{l}\mathrm{u}\left({\mathrm{E}}_{\mathrm{x},0}\right)=6.5\times {10}^{-7} \mathrm{V}/\mathrm{m}, \mathrm{u}\left({\mathrm{E}}_{\mathrm{y},0}\right)=6.5\times {10}^{-7} \mathrm{V}/\mathrm{m}, \mathrm{u}\left({\mathrm{E}}_{\mathrm{z},0}\right)=5.7\times {10}^{-6} \mathrm{V}/\mathrm{m},\\ \mathrm{u}\left({\mathsf{\phi }}_{\mathrm{x},0}\right)=1.7\times {10}^{-3} \mathrm{rad}, \mathrm{u}\left({\mathsf{\phi }}_{\mathrm{y},0}\right)=2.8\times {10}^{-3} \mathrm{rad},\end{array}$$

(22)

for a broad range of angular parameters. The uncertainty of the polarization orientation may be conservatively estimated by assuming that the total electric field uncertainty is perpendicular to the electric field direction. Summing quadratically the uncertainties of the amplitudes of the Cartesian components from Equation (22) and dividing to the amplitude yields an angular uncertainty of 2.1 × 10⁻⁴ rad.

The THz electric field parameters cannot be determined for the following sets of angular parameters $\left(\mathsf{\varphi }={10}^{-11}, \mathsf{\theta }\right), \left(\mathsf{\varphi }=\mathsf{\pi }/2+{10}^{-11}, \mathsf{\theta }\right), \left(\mathsf{\varphi }, \mathsf{\theta }={10}^{-11}\right), \left(\mathsf{\varphi },\mathsf{\theta }=\mathsf{\pi }/2+{10}^{-11}\right)$ because of exceedingly high uncertainties. In this case, the THz electric field may be properly characterized by using a different set of orientations $\left({\mathsf{\alpha }}_i,{\mathsf{\beta }}_i\right)$ of the magnetic field. The lowest uncertainties of the parameters are reached for different orientations of the polarization of the THz wave: $\mathrm{u}\left({\mathrm{E}}_{\mathrm{x},0}\right)=1.4\times {10}^{-8} \mathrm{V}/\mathrm{m}$ for $\left(\mathsf{\varphi }\cong \mathsf{\pi }/30,\mathsf{\theta }\cong 7\mathsf{\pi }/15\right) ,$ $\mathrm{u}\left({\mathrm{E}}_{\mathrm{y},0}\right)=1.4\times {10}^{-8} \mathrm{V}/\mathrm{m}$ for $\left(\mathsf{\varphi }\cong 7\mathsf{\pi }/15,\mathsf{\theta }\cong 7\mathsf{\pi }/15\right) ,$ $\mathrm{u}\left({\mathrm{E}}_{\mathrm{z},0}\right)=1.8\times {10}^{-6} \mathrm{V}/\mathrm{m}$ for $\left(\mathsf{\varphi }\cong \mathsf{\pi }/30,\mathsf{\theta }\cong \mathsf{\pi }/30\right) ,$ $\mathrm{u}\left({\mathsf{\phi }}_{x,0}\right)=6.7\times {10}^{-5} \mathrm{rad}$ for $\left(\mathsf{\varphi }\cong \mathsf{\pi }/30,\mathsf{\theta }\cong \mathsf{\pi }/4\right) ,$ and $\mathrm{u}\left({\mathsf{\phi }}_{y,0}\right)=1.2\times {10}^{-4} \mathrm{rad}$ for $\left(\mathsf{\varphi }\cong 7\mathsf{\pi }/15,\mathsf{\theta }\cong \mathsf{\pi }/4\right) .$ That points at the importance of the choice of the orientation of the magnetic field respective to the polarization of the THz wave and to the potential accuracy with this method. Moreover, even if the uncertainties of the differential standard dynamic polarizabilities are assumed to be zero, the results plotted in Figure 4 remain nearly the same. Therefore, for these measurements, the uncertainty of the electric field characterization is dominated by the experimental uncertainties of the lightshift measurements.

The sensitivity of this method is estimated with the lowest amplitude of the THz wave for which the amplitudes and phases of the Cartesian components of the THz electric field can be determined with uncertainties smaller than the corresponding nominal values. The limit of vector detection of linearly polarized THz waves is estimated here at 350 µV/m, for which the parameters of the THz electric field can be characterized with uncertainties at the $\mathrm{u}\left({\mathrm{E}}_{\mathrm{x},0}{,\mathrm{E}}_{\mathrm{y},0}{,\mathrm{E}}_{\mathrm{z},0},{\mathsf{\phi }}_{\mathrm{x},0},{\mathsf{\phi }}_{\mathrm{y},0}\right)=\left({10}^{-6}{ \mathrm{V}/\mathrm{m}, 10}^{-6}{ \mathrm{V}/\mathrm{m}, 10}^{-4} \mathrm{V}/\mathrm{m}, 1 \mathrm{rad}, 1 \mathrm{rad}\right)$ levels or better. The requirement for the uncertainties of the amplitudes is satisfied alone for a broad range of angular parameters. The requirement for the determination of the phases is satisfied only for narrow ranges of orientation of the THz wave polarization.

4. Discussion and Conclusions

This article proposes a new method to measure a THz electric field by using precision measurements of the lightshifts induced on Zeeman subcomponents of the two-photon rovibrational transition (v, L) = (0, 0)→(2, 0) of cold trapped HD⁺ ions detected by resonance-enhanced multiphoton dissociation. The fractional frequency uncertainty of the infrared transition measurements is estimated at the 10⁻¹² level. An electric field oscillating by 1.3 THz is off-resonantly coupled to the electric-dipole allowed (v, L) = (0, 0)→(0, 1) rotational transition of HD⁺. The lightshifts of the magnetic subcomponents of the ground rotational level may be measured with uncertainties estimated with the molecular ion quantum projection noise limit at 3.5 Hz, using state-selective two-photon rovibrational spectroscopy. The THz electric field is calibrated by comparing the frequency measurements of the lightshifts against a frequency standard and the ab initio predictions of the molecular theory for the energy levels [39,41,42] and for the dipole moments [40]. This approach takes also into account the uncertainties from the theoretical calculations and ensures traceability to the SI units and to a set of fundamental constants, namely the Rydberg constant, the fine structure constant, the proton, deuteron, and electron masses, the proton and deuteron radii, the electric charge, the Planck constant [4], and the deuteron electric quadrupole moment [51].

The lightshift measurement allows scalar sensing of a THz wave. The extension toward polarization measurement is allowed by introducing a directional reference for the molecular ions—that is the quantization axis, which is defined with a static magnetic field. The SI traceability of the magnitude and of the orientation of the magnetic field is allowed by comparing two-photon rovibrational frequencies from Zeeman spectroscopy referenced to a frequency standard, with frequency measurement uncertainty estimated with the molecular ion quantum projection noise limit at 2.57 Hz, with the theoretical predictions for the HD⁺ energy levels. The comparison of theory and spectroscopy results based on the (v, L, F, S, J, J_z)=(0, 0, 1, 2, 2, −2)→(2, 2, 1, 2, 4, 0) transition allows the calibration of the magnetic field magnitude with an accuracy at the 10⁻⁷ T level, and it also allows defining the quantization axis with an angular accuracy at the mrad level.

The electric field of a THz wave can be fully characterized using the response of the HD⁺ molecular ions with respect to the quantization axis direction. The lightshift, depending on the standard components of the THz electric field, is quantified with the dynamic polarizabilities of the HD⁺ energy levels that are calculated here with their uncertainties. This work proposes an algorithm to calculate the values and the uncertainties of the Cartesian components of the THz electric field in the laboratory frame and their relative phases from a set of six independent lightshift measurements performed for two different orientations of the quantization axis. Although similar to the approaches demonstrated with Rabi rate measurements using coherent population transfer [52,53], this method based on lightshift measurements can allow direct referencing of the response to a frequency standard.

The application of this method to a linearly polarized THz wave with an intensity of 1 W/m² (THz electric field amplitude 27.42 mV/m), which is detuned by 1.6 MHz from the (v, L, F, S, J) = (0, 0, 1, 2, 2)→(0, 1, 1, 2, 3) hyperfine component of the fundamental rotational transition of HD⁺, can allow SI calibration of the Cartesian components of the THz electric field with µV/m accuracy and of their phases with mrad accuracy for a broad range of orientations of the THz wave polarization. The accuracy of the phase measurement is two hundred times better than the previous result in microwave electric field phase measurement obtained with the Rydberg atom mixer [21]. The orientation of the polarization is determined with 0.21 mrad accuracy, which is forty times better than the vector microwave electrometry result obtained with Rydberg atom spectroscopy [20]. The lowest THz electric field vector that can be completely characterized with this method has an amplitude estimated at 350 µV/m. The detection limit in scalar THz electrometry at 11 µV/m estimated for HD⁺ ion spectroscopy [54] represents a 16 times improvement to the result obtained in Rydberg atom-based scalar microwave electrometry [18]. That points to the potential of lightshift measurements for THz electric field sensing—a method extensible to a number of cold trapped molecular ion species, to the optical fields, and to the magnetic components, which is compatible with the quantum enhancement techniques [32,33]. The significant improvements in performances allowed by this method compared to the traditional techniques can be of general relevance to the field of the electric field metrology.