text:
\section{Introduction}

Over the next couple of years, a new generation of wide-field optical galaxy surveys will be coming online, starting with the ESA Euclid mission \citep{Laureijs2011} and the Vera C. Rubin Observatory Legacy Survey of Space and Time \citep{Ivezic2019}, soon followed by the Roman Space Telescope \citep{Spergel2015}. At the core of the cosmology science that these surveys will enable is the measurement of the weak gravitational lensing effect, a very subtle shearing of the apparent images of distant galaxies in the presence of massive structures along the line of sight. On cosmological scales, this effect, also known as cosmic shear, is an invaluable probe of the total matter content of the Universe, and will play a key role in further constraining and testing our cosmological model over the next decade.

Despite over 20 years of research in this field, unbiased measurements of gravitational shear from observed galaxy images remains one of the most challenging aspects of the scientific exploitation of these surveys. The prevailing strategy for measuring shear from data has been to build estimators of the ellipticity of galaxy light profiles, due to the fact that one can in theory relate how a galaxy ellipticity $e$ changes under a gravitational shear $\gamma$, as $e_\text{obs} \simeq e_{i} + \gamma$. Under the assumption that galaxies should have no preferential intrinsic ellipticity, a simple estimator for $\gamma$ can therefore be constructed from $\hat{\gamma} = \langle e_\text{obs} \rangle = \langle e_i \rangle + \gamma$. However, to this day, no shear measurement algorithm has proven to be unbiased.

In fact, state-of-the-art approaches for shear measurement \cite{huff2017metacalibration, sheldon2017practical} no longer attempt to build an unbiased shear estimator from images, but rather calibrate the shape measurement methods estimating a multiplicative and an additive bias. There are many factors contributing to these biases, but a fundamental issue comes from the fact that the ellipticity $e$ is not a well defined quantity for arbitrary galaxy light profiles.

In this paper we propose a new paradigm for cosmic shear estimation where we no longer rely on ellipticities anymore. We propose to forward model the galaxy images at the pixel level, similarly to \cite{schneider2015hierarchical}, but using non-parametric galaxy light profiles. The main idea is to train a deep generative model to learn galaxy morophologies, which can then be used as part of a hierarchical Bayesian model, a hybrid deep learning and physical model, simulating the observations and used to infer the shear.

The purpose of this paper is first to demonstrate that using parametric light profiles as part of the forward model induces a bias in the shear estimation, which we call a \textit{model bias}, and then to show that this bias can be removed by using a deep generative model for the morphology. We show in this preliminary study, that using a hybrid deep learning-physical Bayesian forward model can not only solve the model bias, but also provide uncertainty quantification to the shear estimation.

\begin{figure}[h]
    \centering
    \begin{tabular}{cc}

        \includegraphics[width=0.47\textwidth]{assests/sersicfit1.pdf} & \includegraphics[width=0.47\textwidth]{assests/deepgal1.pdf} \\

        (a) \small{Realistic galaxies fitted by a parametric profile} & (b) \small{Realistic galaxies fitted by a deep morphology model}
    \end{tabular}
    \caption{The first row of images correspond to observed galaxies. The second row corresponds to the light profiles Maximum a Posteriori fits, convolved with their PSFs. The third row shows  residual images between the observations and the model fits. Column (a) is using Sérsic model as light profile and column (b) is using \texttt{deepgal} model as light profile.
    }
    \label{fig:MAP}
\end{figure}

\section{Forward modeling weak lensing observations at the pixel level} \label{sec:models}

From a Bayesian perspective, making the joint inference of cosmic shear and galaxy morphologies means estimating the \textit{joint posterior} distribution $p(\mathcal{G}, \gamma | \mathcal{D})$, where $\mathcal{G} = \{z_i\}_{i=1\dots N}$ is the set of morphology parameters and $\mathcal{D} = \{d_i\}_{i=1\dots N}$ is the set of observed data. The Bayes identity tells us:
\begin{equation}\label{eq:bayes}
    p(\mathcal{G}, \gamma | \mathcal{D}) \propto p(\mathcal{D} | \mathcal{G}, \gamma) p(\mathcal{G}, \gamma),
\end{equation}
where $p(\mathcal{D} | \mathcal{G}, \gamma)$ is the data \textit{likelihood} given by the forward model described bellow. Since the intrinsic morphology of galaxies are independent of the cosmic shear, then $p(\mathcal{G}, \gamma) = p(\mathcal{G})p(\gamma)$, where $p(\mathcal{G})$ and  $p(\gamma)$ are respectively the \textit{prior} distributions of the morphological parameters and of the cosmic shear.

We show the probabilistic graphical model in \autoref{fig:pgm}. Nodes denote random variables and arrows conditional relation between them. White nodes corresponds to random variable on which we do the inference on, conditioned on the shaded nodes corresponding to observed variables. Dot nodes correspond to fixed variables, e.g. neural network weights $\theta$. In this paper we compare two forward models: (1) using Sérsic law  and (2) using the \texttt{deepgal} model to generate the galaxy light profile. Once the galaxy light profiles are drawn, we shear images with a constant value $\gamma$ shared by all the galaxies. Each sheared galaxy $i$ is then convolved by its Point Spread Function (PSF) $\Pi_i$. Assuming Gaussian noise corruption, the resulting image is therefore evaluated by a Gaussian likelhood of variance $\sigma_\text{pix}^2$.

\begin{wrapfigure}[12]{r}{0.4\textwidth}
    \vspace{-25pt}
    \centering
    \includegraphics[width=0.35\textwidth]{assests/pgm_deepgal.pdf}
    \vspace{-10pt}
    \caption{Probabilistic Graphical Model.
    }
    \label{fig:pgm}
\end{wrapfigure}
\textbf{Sérsic model~~}
The light profile of the Sérsic model is given by the following parametrization:
$I(r) = F \times I_0 \exp ( -b_n [( \frac{r}{r_\text{hlr}})^{\frac{1}{n}} -1] ),$
where $r_\text{hlr}$ is the half-light-radius, $I_0$ is a normalization factor such that the integration of $I$ over all the space is $F$, which corresponds to the flux of the galaxy. $n$ is called the Sérsic index and $b_n$ is a dimensionless constant depending on $n$. The galaxy is then sheared with parameters $e=(e_1, e_2)$ corresponding to its ellipticity. Finally, since the real galaxy may not be perfectly centered in the postage stamp, we add two shift parameters $s_x$ and $s_y$. Therefore, each galaxy $i$ has a set of parameter $z_i = \{n, r_\text{hlr}, F, e_1, e_2, s_x, s_y\}_i$. Prior distributions and parameters follow empirical distributions of the COSMOS shape catalog and are described in \autoref{sec:cosmos-prior}.

\textbf{\texttt{Deepgal} model~~} \texttt{deepgal} is a VAE-based model, proposed in \citep{lanusse2021deep}, generating complex galaxy light profiles. In this model, the observations are generated following a random process mapping a latent space representation $z_i \in \mathbb{R}^{16}$ to the observed image $d_i$ following the probabilistic model illustrated in \autoref{fig:pgm}. The neural network weights $\theta$ of the VAE are trained maximizing the Evidence Lower Bound (ELBO) of the data. Moreover, another generative model was trained on top of the VAE to learn latent space distribution using Masked Auto Regressive Flows \cite{papamakarios2017masked}. This model was trained on the HST/ACS COSMOS survey postage stamps.

In the weak lensing regime, the shear is very small, therefore we can assume that it does not change the morphological statistics of the training sample. Thus, the model can be trained once, a priori, over a large catalog of galaxies and be used to perform the joint inference of the latent variables and the cosmic shear. For both light profile models, we denote the generated image $\mathcal{I}_0$.

\textbf{Cosmic shear~~} In this model every galaxy image is sheared by a constant value $\gamma = \gamma_1 + i \gamma_2$ shared by all the galaxies. The image shearing transform, denoted $\mathcal{I}_0 \otimes \gamma$, is given by the following pixel coordinate change: $\big(\begin{smallmatrix}
  x^\prime \\
  y^\prime
\end{smallmatrix}\big) = S(\gamma) \big(\begin{smallmatrix}
  x \\
  y
\end{smallmatrix}\big)$,
such that $S(\gamma) = \frac{1}{\sqrt{1-|\gamma|^2}}
                    \big(\begin{smallmatrix}
  1+\gamma_1 & \gamma_2 \\
  \gamma_2 & 1-\gamma_1
\end{smallmatrix}\big)$.

\textbf{Convolution with the PSF~~}
The Point Spread Function (PSF) model is assumed to be known and exact for each galaxy. We draw the PSF evaluated at the galaxy $i$ position, denoted by $\Pi_i$, in Fourier space and then multiply it to the Fourier image of the light profile before returning in the real space, resulting in an image $\mathcal{I}$, such that $\mathcal{\tilde I} =  (\mathcal{\tilde I}_0 \otimes \gamma) \cdot \tilde \Pi$.

\textbf{Noise corruption~~}We assume a known fixed gaussian noise identical for all the galaxies. For simplicity, all of the galaxies in our observations are corrupted by a constant Gaussian noise of variance $\sigma_\text{pix}^2$, as well as for our forward model. In a more realistic scenario where we would have information about the noise correlation properties for each postage stamp, one could consider a noise model per galaxy. The likelihood is therefore given by

\begin{equation}
    p(\mathcal{D} | \mathcal{G}, \gamma) \propto \prod_{i=1,\dots, N} \mathcal{N}\left( d_i \Big|  \mathcal{F}^\dagger\left[ \left(\mathcal{\tilde I}_{0,i} \otimes \gamma \right) \cdot \tilde \Pi_i \right], \sigma^2\right), ~~~\mathcal{F}^\dagger~ \text{is the inverse Fourier operator}.
\end{equation}

\vspace{-10pt}

\section{Illustration on realistic data}

\begin{figure}[h]
    \centering

    \includegraphics[width=\textwidth]{assests/biases5.pdf}
    \caption{Shear marginal posterior contours against ground truth $\gamma=(-0.03, 0.03)$. The left panel shows the posterior using the Sérsic light profile over the parametric observations. The middle panel shows the posterior using the Sérsic light profile over the realistic galaxies. The right panel shows the posterior using the \texttt{deepgal} model over the realistic galaxies.
    }
    \label{fig:bias}
\end{figure}

\vspace{-10pt}

\subsection{Simulating observations} \label{sec:data}

We designed a toy model where galaxies are assumed to be already detected, deblended and centered in postage stamps, so that we only need to simulate the light profile, shearing transform and convolution with the PSF. In order to demonstrate the model bias, we distinguish two sets of observations. (1) generated with parametric galaxy light profiles using a Sérsic model and (2) using realistic galaxies based on the galaxy postage stamps from the COSMOS HST Advanced Camera for Surveys (ACS) field \citep{koekemoer2007cosmos, scoville2007cosmic_a, scoville2007cosmos_b}. We used the postage stamps from \cite{leauthaud2007weak} and Sérsic fit parameters compiled in \cite{mandelbaum_rachel_2012_3242143}. Because the unconvolved galaxies are very noisy, we used the Auto-Encoder from \citep{lanusse2021deep} to get denoised version of unconvolved galaxies, therefore this is not a demonstration on real galaxies yet. However, the purpose of this paper is to demonstrate the model bias, so we find that realistic galaxy morphologies that cannot be captured by analytic profiles are sufficient to illustrate the argument.

Each set of observation is a field of 256 galaxies, generated from the COSMOS catalog, with the following cuts: $n > 0.4$, $r_\text{hlr} > 0.2$ and $23.5 < i^+ < 25.2$ so rather faint galaxies. That way the entire galaxy is in the postage stamp, is not too peaky and not too resolved to be close to a Sérsic model. Every postage stamps were corrupted by a Gaussian noise of standard deviation $\sigma_\text{pix}=0.01$. In order to simulate shearing, we removed any shear present in the selected sample from the COSMOS catalog, by adopting the same strategy than in the GREAT3 Challenge \citep{mandelbaum2014third}, i.e. doubling each galaxy, and rotating the doubled galaxies by 90°, as well as their PSF, to cancel the overall sample shear and then apply the same constant shear to every galaxy $\gamma=(-0.03, 0.03)$.

\subsection{Approximate inference with mean-field VI}

In order to estimate the posterior distribution $p(\mathcal{G}, \gamma | \mathcal{D})$, we relied on approximate inference using mean-field Variational Inference (VI), where we approximate the posterior distribution using a parametric distribution $q_\phi(\mathcal{G}, \gamma) = q_{\phi_\gamma}(\gamma) \prod_i\prod_j q_{\phi_{i,j}}(z_{i,j})$,
parameterized by $\phi$. We used Normal distributions $q_{\phi_{i,j}}$ to approximate the posteriror of parameters $i$ of the $j^\text{th}$ galaxy, so the parameters $\phi_{i,j} = (\mu_{i,j}, \sigma_{i,j})$ correspond to means and standard deviations. Similarly, $q_{\phi_\gamma}$ is a 2d-diagonal multivariate Gaussian. The parameters $\phi_{i,j}$ and $\phi_\gamma$ are then optimized maximizing the ELBO:
\begin{equation}
    \phi^* = \arg \max_\phi ~~ \mathbb{E}_{(\mathcal{G}, \gamma) \sim q_\phi} \left[ \log p(\mathcal{G}, \gamma, \mathcal{D}) - \log q_\phi(\mathcal{G}, \gamma) \right]
\end{equation}
This remains a very simplistic approach, which does not assume any correlation between variables. A more accurate model would condition the morphological parameters on the shear. In principle, the most accurate method would be to run a Markov Chain Monte Carlo, but several experiments showed that sampling based-methods are too computationally expensive for such high-dimensional problem.
\subsection{Identifying and solving the model bias}

We ran three sets of experiments to demonstrate the model bias. First fitting the fully parametric model to the parametric observations. Then fitting the same parametric model to the realistic galaxies and then fitting the hybrid model to the realistic galaxies. We used a Sérsic law as light profile as proposed in \cite{schneider2015hierarchical}. One could argue that there exist more complex parametric profiles such as the bulge+disk model which could reduce the aforementioned bias, but we propose here a model which does not rely on ellipticities and in principle can capture any complexity of galaxy morphologies.

First, as illustrated in \autoref{fig:bias} (a), the inferred shear using the parametric model on the parametric observations is unbiased. This was indeed expected since the observations and the model are exactly the same. Secondly, \autoref{fig:bias} (b) shows that the parametric model induces a bias in the shear estimation when the galaxy morphologies are no longer elliptical. Besides, one can observe the non-zero residuals in the \autoref{fig:MAP} fit to real galaxies, illustrating the model discrepancy when applied to complex morphologies. Finally, \autoref{fig:bias} (c) shows that when we substitute the parametric light profile model by the \texttt{deepgal} model, the shear estimation becomes unbiased. One can notice that the uncertainties using the generative model are smaller than for the Sérsic model, but these contours are not representative of the error bars due to the mean-field approximation.

\section{Discussion and conclusion}

We have presented in this work the necessity of an accurate light profile model for forward modeling of the cosmic shear. We have shown that a hybrid physical and deep learning HBM has the potential to solve the model bias, being able to capture the complex morphology of galaxies and simulating physical processes such as shearing transform and convolution by the PSF. We identify many directions of research towards a real cosmological analysis by forward modeling. On the one hand, many assumptions could be relaxed by building a more complex hierarchical model, e.g. parametrizing the PSF model and jointly inferring it along with the shear. Concerning the detection and deblending of galaxies, similar work have already been proposed, e.g. BLISS \cite{hansen2022scalable}. On the other hand, accurate uncertainty quantification could be enabled by scaling MCMCs for HBMs or using neural density estimator to capture the shape of the joint posterior distribution.

\section{Broader impact}

The methodology described here finds many useful applications within and outside of cosmology. Hierarchical Bayesian Models are very efficient tools for statistical inference but can suffer form model misspecification. Combining deep generative models and physical models seems to be a promising approach to solve this problem. We believe this work does not entail any negative consequences or ethical issues.

\bibliographystyle{unsrt}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{intro}
The quantum magnet \lirefx{Ho} has been shown to exhibit a variety of quantum many-body phenomena, such as quantum annealing \cite{brooke1999}, long-lived coherent oscillations \cite{ghosh2002}, long-range entanglement \cite{ghosh2003}, quantum phase transitions \cite{ronnow2005} and high-$Q$ nonlinear dynamics \cite{silevitch2019}. Importantly, knowledge about the hyperfine (HF) interactions and the lowest-energy crystal-field (CF) states is necessary for the understanding of the material's properties.

In addition, rare-earth doped crystals are promising candidates for quantum information applications \cite{ortu2018,kindem2020,raha2020,grimm2020}. Among the many potential host materials, isotopically-pure \lirefx{Er} exhibits long electronic coherence times~\cite{kukharchyk2018}. To control the electro-nuclear degrees of freedom \eg \ for quantum information processing, precise knowledge of the underlying electronuclear Hamiltonian is required.

Owing to limitations in resolution and other experimental challenges, comprehensive high-resolution measurements of the transitions within all low CF levels of \lirefx{Ho} have not been available so far~\cite{magarino1980,kjaer1989,babkevich2015,matmon2016}, although we have recently combined optical comb synthesis and a software controlled  modulator to obtain ultra-high resolution data for the HF-split lowest CF excitation near 6.8\,$\textrm{cm}^{-1}$~\cite{hermans2020}. Here we integrate theory together with low-temperature terahertz time-domain spectroscopy (TDS) and synchrotron-based ultra-high resolution Fourier transform infrared (FTIR) spectroscopy (Sec.~\ref{method}) to re-examine the transitions between the three lowest-lying CF levels of the $^5$\textit{I}$_8$ ground-state manifold of \lirefx{Ho}.

Motivated by the precisely and unambiguously resolved HF splitting, we expand previous treatments of the HF interaction in CF states (\eg~\cite{elliot1972,matmon2016}) up to second order in the dipolar HF interaction, and to first order in the quadrupolar coupling in Sec.~\ref{CF transitions and HF}. There, we discuss the measured data, in particular the transition energies including isotopic shifts due to $^6\mathrm{Li}$ and $^7\mathrm{Li}$.
We extract CF parameters by combining our data with CF energy measurements of the $^5$\textit{I}$_8$ manifold from Ref.~\cite{matmon2016}, which enables us to refine the dipolar HF interaction constant $A_J$ in Sec.~\ref{sec: CF parameters}. Based on the resulting CF parameters, we predict all $^5$\textit{I}$_8$ CF energies and their magnetic moments. The high instrumental resolution allows us also to determine the quadrupolar HF constant $B$. Using $B$ and $A_J$, we infer non-equidistant HF corrections of the three CF levels involved in our measurements, including the ground state. We provide an approximation of these HF corrections based on our measurement results, which corroborate our numerical simulation. We conclude in Sec.~\ref{sec: CF parameters} by discussing the implications of non-equidistant HF corrections on unambiguous addressing of specific HF transitions \eg \ for quantum information applications. In the appendix~\ref{Sec: Refractive Index} we provide a refractive index measurement of \lirefx{Ho} from $70$ to $5$\,$\textrm{cm}^{-1}$ \ corresponding to $\sim2$~THz to the sub-THz range; these data are useful for planning the design of future optical experiments and devices. A summary and an outlook are found in Sec.~\ref{conclusion}.

\section{Experimental setup}
\label{method}

\subsection{Sample}
For an overview of the physical properties of \lirefx{Ho}, we refer to Refs.~\cite{gingras2011,matmon2016}. We study three commercially-available \lirefx{Ho} single crystals at low doping concentrations of \mbox{$x=1\%$}, $0.1\%$, and $0.01\%$. The crystal dimensions along the light propagation direction are chosen such that transmission is optimized for each $x$. Samples were mounted on the cold finger of a continuous-flow liquid-helium cryostat. The THz light was linearly polarized. The sample was oriented with the crystallographic $c$-axis parallel to the magnetic field component, while the light's propagation direction was perpendicular to $c$. All reported temperatures denote the nominal values on the cryostat cold finger.

\subsection{Experimental methods}
We use two different measurement methods: First, TDS was conducted on \lirefx{Ho} ($x=0.1\%$) for wavenumbers $\tilde{\nu}<10$\,$\textrm{cm}^{-1}$ \ (300\,GHz), as well as for refractive-index measurement of the $x=1\%$ crystal for $\tilde{\nu}\le70$\,$\textrm{cm}^{-1}$ \  (2.1~THz).
Figure~\ref{fig:TDSsetup} shows a schematic of the custom experimental setup, which is based on an 800\,nm laser, delivering 100\,fs pulses at 80\,MHz repetition rate. The beam is split, directing 250\,mW through a variable delay line. This fraction of the laser is focused onto a low-temperature grown \chem{GaAs} photo-conductive emitter with a 100\,$\mu$m electrode gap (biased at 100\,V, 7.3\,kHz) that generates a linearly polarized single-cycle THz pulse. The THz pulse is then collected from the back of the emitter substrate with a Si hyper-hemispheric lens and focused on the sample with a parabolic mirror. Thereafter, the transmitted beam is refocused onto a 2\,mm thick \chem{ZnTe} crystal for electro-optic sampling. In this detection scheme, the THz branch is overlapped with the 800~nm branch. As a function of delay time, the polarization change of the transmitted 800~nm light is then proportional to the instantaneous THz field in the ZnTe crystal. The signal is measured using balanced photo-diodes and a lock-in amplifier referenced to the emitter bias frequency. Fourier transforms of the delay scans then yield the spectra.

\begin{figure}
\includegraphics[width=\columnwidth]{Figures/setupTHzTDS.pdf}
\caption{Schematic of the time-domain spectroscopy setup showing the femtosecond (fs) laser, the beam-splitter (BS), a photo-conductive (PC) antenna, the sample (S), a ZnTe electro-optic (EO) sampling crystal, a quarter-wave plate (QWP), a Wollaston prism (WP) and balanced photo-diodes (BD) for detection.}
\label{fig:TDSsetup}
\end{figure}

Second, ultra-high resolution FTIR spectroscopy was conducted on \lirefx{Ho} ($x=0.01\%$) for \mbox{$\tilde{\nu}>15$\,$\textrm{cm}^{-1}$}~(450\,GHz) using a custom-built Bruker FTIR spectrometer with 0.00077\,$\textrm{cm}^{-1}$ \ (23\,MHz) resolution. A He-flow cryostat for low-temperature measurements was fitted to the spectrometer. The THz source is the high-brilliance, highly collimated far-infrared (FIR) radiation from the Swiss Light Source synchrotron at the Paul Scherrer Institut, Switzerland. Reference~\cite{albert2011b} provides more details about this FTIR setup. The unique combination of a low-temperature, ultra-high resolution spectrometer and FIR synchrotron radiation allowed us to measure the absorbance spectra with a resolution of up to $10^{-3}$\,$\textrm{cm}^{-1}$, which is more than an order of magnitude higher than previously reported ~\cite{karayianis1976,christensen1979,matmon2016}.

The THz response of the holmium ions (\ion{Ho}) in the \lyf \ matrix is characterized by referencing the sample absorption at low temperature to a higher temperature measurement. This ensures that both the background absorption of the crystal host and temperature-independent reflections from the sample and the experimental setup are removed. Therefore, we show absorbance spectra $A(\tilde{\nu})=\mathrm{Log}_{10}[I_0(\tilde{\nu})/I(\tilde{\nu})]$ as a function of wavenumber $\tilde{\nu}$ {[$\textrm{cm}^{-1}$]}, with $I(\tilde{\nu})$ $(I_0(\tilde{\nu}))$ denoting the wavenumber-dependent sample (reference) transmission.

\section{Crystal field transitions with hyperfine interactions}
\label{CF transitions and HF}

Our high-resolution setups enable us to resolve the HF structure of the measured CF states to high precision; analysis methods which take advantage of this structure are described in Ref.~\cite{beckert2020}. We turn now to the theoretical understanding of the HF corrections to the measured CF states. In this paper we denote a transition from an initial CF state $i$ to a final state $f$ by $i\rightarrow f$. Further, we label the $^5$\textit{I}$_8$ ground state manifold states $8.n$ according to their CF energy $E_n$: the ground state ($8.1,E_1=0$) is a doublet (under time-reversal symmetry) and carries $\Gamma_{3,4}$ symmetry, the first excited (8.2) and second excited (8.3) states have $\Gamma_{2}$ symmetry at $E_2=6.8$\,$\textrm{cm}^{-1}$ \ and $E_3=23.3$\,$\textrm{cm}^{-1}$, respectively. We denote the CF symmetries (irreducible representations) by $\Gamma_j,~j\in\{1,2,3,4\}$, using standard conventions.
Individual HF states are labelled as $\ket{8.n^\sigma, m_z} \equiv \ket{8.n^\sigma} \otimes \ket{m_z}$, where $\sigma = -1$ ($\sigma=+1$) denotes the $\Gamma_3$ ($\Gamma_4$) state if the $n$-th level belongs to a doublet. $m_z$ is the nuclear spin projection onto the crystallographic $c$-axis.

\subsection{HF interaction in perturbation theory}
Within the lowest $J$-multiplet, the electrons of each \ion{Ho} ion ($J=8$) couple to their nuclear spin ($I=7/2$) via the dipolar and quadrupolar HF interactions
\begin{equation}
\begin{split}
    H_{\mathrm{HF}} =& H_{\mathrm{HF, dip}} + H_{\mathrm{HF, quad}} \\
    =& A_J \, \vec{J} \cdot \vec{I} + \frac{B}{2I(2I-1)J(2J-1)} \bigg( 3(\vec{J} \cdot \vec{I})^2 \\
    &+ \frac{3}{2} (\vec{J} \cdot \vec{I}) - I(I+1)J(J+1) \bigg),
\end{split}
\label{eq: HF interaction}
\end{equation}
with the dipolar and quadrupolar coupling constants $A_J$ and $B$, respectively. Below we consider the effects of $A_J$ up to second order and $B$ to first order because of the relative size of these terms. We neglect HF corrections due to coupling of the nuclear electric quadrupole moment to the electric field gradient. Using the literature value in Ref.~\cite{popova2000}, this effect is estimated to be an order of magnitude smaller than the terms in Hamiltonian~(\ref{eq: HF interaction}).

Using perturbation theory, the ground-state energy corrections $\delta_{8.1^\sigma, m_z}$ of the states $\ket{8.1^\sigma, m_z}$ are
\begin{equation}
\begin{split}
    &\delta_{8.1^+,+m_z} = \delta_{8.1^-,-m_z}= A \braket{8.1^+|J_z|8.1^+} m_z\\
    &+ \sum_{ j\in \Gamma_{1}}  \frac{A_J^2}{4 \Delta E_{1j}} \bigg[  |\braket{8.j|J_+|8.1^+}|^2 \big( I(I+1)-m_z(m_z+1) \big) \bigg] \\
    &+ \sum_{j \in \Gamma_{2}}  \frac{A_J^2}{4 \Delta E_{1j}} \bigg[  |\braket{8.j|J_-|8.1^+}|^2  \big( I(I+1)-m_z(m_z+1) \big) \bigg]\\
    &+ \sum_{\substack{j \in \Gamma_{3,4}\\ j \neq 1}} \frac{A_J^2}{\Delta E_{1j}} \bigg[ |\braket{8.j^+|J_z|8.1^+}|^2 m_z^2  \bigg]\\
    &+ \frac{B \braket{8.1^+|3J_z^2-J(J+1)|8.1^+}}{4I(2I-1)J(2J-1)}  (3m_z^2-I(I+1)) ,
\end{split}
\label{eq: HF corr doublet}
\end{equation}
and the corrections of the first two excited electronic states ($n=2,3$) are
\begin{equation}
\begin{split}
    &\delta_{8.n, \pm m_z} =  \sum_{\substack{j \in \Gamma_{2}\\ j\neq n}}  \frac{A_J^2}{\Delta E_{nj}}  |\braket{8.j|J_z|8.n}|^2 m_z^2 \\
    &+ \sum_{j \in \Gamma_{3,4}} \frac{A_J^2}{2 \Delta E_{nj}}  \bigg[ |\braket{8.j^+|J_+|8.n}|^2  \big( I(I+1)-m_z(m_z+1) \big) \bigg]\\
    &+ \frac{B \braket{8.n|3J_z^2-J(J+1)|8.n}}{4I(2I-1)J(2J-1)}  (3m_z^2-I(I+1)) .
\end{split}
\label{eq: HF corr singlet}
\end{equation}
Here $\Delta E_{nj}=E_{n}-E_{j}$ is the energy difference between the CF levels $\ket{8.n}$ and $\ket{8.j}$. The sums run over all CF states $\ket{8.j}$ carrying the irreducible representations $\Gamma_i,~i\in\{1,2,3,4\}$. From now on, we use the abbreviation $\lambda_n/2$ for the prefactor of the HF corrections $\propto m_z^2$ of the CF states $8.n$.

These perturbative corrections are sufficient to interpret the HF spectrum of the 8.1, 8.2 and 8.3 states, where $m_z$ denotes the projection of the nuclear spin in the unperturbed electron-nuclear wavefunction. Due to the absence of external magnetic fields, then by Kramers' theorem all HF states are doubly degenerate with their time-reversed state (under time-reversal: $m_z \to -m_z$, $\sigma \to -\sigma$). The electronic doublet 8.1, which is Ising-like with a moment along the crystallographic $c$-axis (due to the $S_4$ site symmetry~\cite{gingras2011}), experiences a dominant first order shift $\propto A_J m_z$ that leads to an equidistant HF splitting into eight HF Kramers doublets. In the lowest ($m_z=-7/2$) and highest ($m_z=+7/2$) of these HF states the electronic and magnetic moments are anti-aligned and aligned, respectively. The singlets do not undergo a first-order HF shift in $A_J$ due to their vanishing moment. Within a single CF state, the equidistance of the HF energies is broken by the second-order terms in $A_J$ and first-order term in $B$, all leading to corrections $\propto m_z^2$. These corrections determine the relative order of the $m_z$ states within a singlet. For the states 8.2 and 8.3, the relative order is reversed. The dominant correction due to the small energy denominator in Eq.~(\ref{eq: HF corr singlet}) comes from the mutual repulsion of the CF states caused by the dipolar HF interaction. An illustration of the HF levels of the 8.2 and 8.3 states is shown in Fig.~\ref{FigTransition82to83}.

\subsection{Experiments}
We measure the transmission of the $8.1\rightarrow8.2$ magnetic dipole transition in \lirefx{Ho} ($x=0.1\%$) at a temperature of $T=2.9$\,K by TDS with an instrument resolution of 0.017\,$\textrm{cm}^{-1}$ \ (500\,MHz). The absorbance is shown in Fig.~\ref{specplots}(a), where we directly resolve an eight-fold, approximately equidistant HF splitting of $\sim0.146$\,$\textrm{cm}^{-1}$ \, which reflects the dominant linear HF shift of the ground state doublet 8.1.
The deviation of the individual line intensities from a Boltzmann distribution (\cf \ Refs.~\cite{babkevich2015,matmon2016}) originates from sample- and setup-specific systematic errors such as residual interference of optical components. The extracted Gaussian full-width-at-half-maximum (FWHM) of a single HF line is $0.017\pm0.001$\,$\textrm{cm}^{-1}$ \ and thus instrument-resolution limited.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{Figures/spectra.pdf}
\caption{Absorbance spectra of \lirefx{Ho} of (a) the $8.1\rightarrow8.2$ (TDS, $x=0.1\%$, $T=3$\,K) and (b) the $8.1\rightarrow8.3$ (FTIR, $x=0.01\%$, $T=3.5$\,K) transitions with conserved $m_z$.}
\label{specplots}
\end{figure}

The absorbance of the $8.1\rightarrow8.3$ magnetic dipole transition of \lirefx{Ho} \mbox{($x=0.01\%$)} was measured at $T=3.5$\,K with FTIR spectroscopy and 0.001\,$\textrm{cm}^{-1}$ (30\,MHz) resolution. The absorbance spectrum is shown in Fig.~\ref{specplots}(b), also revealing the eight-fold CF level splitting. The HF lines are nearly equidistant with a spacing of $\sim0.146$\,$\textrm{cm}^{-1}$. The ultra-high resolution of the FTIR spectrometer allows for a closer inspection of a single HF line. Figure~\ref{detspec} shows the sixth HF peak at $23.527\pm0.001$\,$\textrm{cm}^{-1}$ \ in more detail. An asymmetry towards larger wavenumbers is apparent, and is best explained by the isotopic splitting effect due to the natural abundance of $^6$Li ($7.6 \%$) and $^7$Li ($92.4 \%$). A finite number $\iota$ of lighter $^6$Li atoms which substitute the more abundant $^7$Li in the immediate neighborhood of a \ion{Ho} ion was shown to lead to a slight shift in the crystal field parameters of \ion{Ho}, attributed to a difference in the zero point motion of $^6$Li and $^7$Li~\cite{agladze1991,shakurov2005}. This leads to additional peaks in the absorbance spectrum from Ho$^{3+}$ ions with different numbers of less abundant $^6$Li neighbors. Their intensities decrease exponentially with $\iota$~\cite{agladze1991,shakurov2005}, reflecting the Bernoulli distribution of the number of $^6$Li neighbors. We only take the two strongest peaks $\iota=\{0,1\}$ into account, as peaks corresponding to $\iota>1$ were not observed. By fitting two Gaussians, we find an isotopic splitting of $0.0098\pm0.0004$\,$\textrm{cm}^{-1}$ \ and a Gaussian FWHM of $0.0090\pm0.0001$\,$\textrm{cm}^{-1}$ \ for the individual peaks with the errors extracted from the covariance matrix. These findings are in agreement with the previously reported values of $0.0105\pm0.0015$\,$\textrm{cm}^{-1}$~\cite{shakurov2005}. HF line energies of the $8.1\rightarrow8.3$ transition are always referred to the center of the dominant $\iota=0$ peak.

\begin{figure}[ht]
\centering
\includegraphics[width=0.6\columnwidth]{Figures/isoplot.pdf}
\caption{Close-up of the asymmetric sixth $8.1\rightarrow8.3$ HF peak of \lirefx{Ho} \mbox{($x=0.01\%$)} at $T=3.5$\,K. The green line is a fit to the main peak (blue), attributed to the majority of Ho ions having $\iota=0$ $^6$Li neighbors, and the smaller second peak (red) that is shifted by isotopic splitting due to $\iota=1$ $^6$Li neighbors. A sinusoidal background owing to interference effects was subtracted w.r.t. the data shown in Fig.~\ref{specplots}.}
\label{detspec}
\end{figure}

Next we completed FTIR measurements of the $8.2\rightarrow8.3$ magnetic dipole transition of \lirefx{Ho} \mbox{($x=0.01\%$)} with 0.002\,$\textrm{cm}^{-1}$ \ resolution. The temperature was set to $T=9$\,K, to thermally populate the 8.2 state. The inset of Fig.~\ref{FigTransition82to83} shows the respective absorbance with a Lorentzian fit. The HF corrections $\propto m_z^2$ in Eq.~(\ref{eq: HF corr singlet}) lead to an observable difference in the transition energies of the individual $m_z$ states. The HF levels are also illustrated (not to scale) in Fig.~\ref{FigTransition82to83}. We fit the absorbance spectrum of the $8.2\rightarrow8.3$ transition with four Lorentzian profiles, taking the degeneracy of $\pm m_z$ into account. We allowed for different intensities and peak frequencies, but imposed an identical linewidth, which we found to be $0.013\pm0.001$\,$\textrm{cm}^{-1}$.
Beyond a 0.008\,$\textrm{cm}^{-1}$ \ constant offset, we obtain results that are consistent with the difference measured at the $8.1\rightarrow8.3$ and $8.1\rightarrow8.2$ transitions. We attribute the offset partly to the lower resolution of the TDS setup (0.017\,$\textrm{cm}^{-1}$) and systematic differences between the two experimental setups.

\begin{figure}[t]
\centering
\includegraphics[width=1\columnwidth]{Figures/quadratic16fitplot.pdf}
\caption{\noindent HF shifts of the $8.2\rightarrow8.3$ transition as a function of the nuclear spin  $m_z$. The shifts depend quadratically on $m_z$, as evidenced by the fit (orange). Note the degeneracy of $\pm m_z$ (see Table~\ref{etab}). The top left inset shows the $8.2\rightarrow8.3$ absorbance raw data of \lirefx{Ho} ($x=0.01\%$) at $T=9$\,K (blue) and a fit to Lorentzian profiles (red). The vertical lines denote the peak center positions and correspond to the data in the main figure. Also shown on the bottom right is the energy level diagram of the four observable $8.2\rightarrow8.3$ HF-split CF transitions.}
\label{FigTransition82to83}
\end{figure}

We summarize the HF transition energies in Table~\ref{etab}. Note that the $8.1\rightarrow8.3$ ($x=0.1\%$) transitions, obtained with FTIR, exhibit smaller uncertainties than the $8.1\rightarrow8.2$ ($x=0.01 \%$) transitions, since the latter was measured with lower instrument resolution of the TDS setup. Owing to the significant line overlap of the $8.2\rightarrow8.3$ ($x=0.01 \%$) transition data, the respective uncertainties extracted from the fit covariance matrix amount to $\le 0.003$\,$\textrm{cm}^{-1}$. Reference~\cite{koenz2003} shows that increasing the rare-earth concentrations up to $10\%$ does not noticeably affect the CF energies, which justifies a direct comparison of the $x=0.01\%$ and $0.1\%$ CF energies.

\begin{table}[t]
\begin{tabular}{@{}c|l|ccc@{}}
\toprule
HF Index~~ &~~$m_z$~~ & ~~$8.1^{+}\rightarrow8.2$~~ & ~~$8.1^{+}\rightarrow8.3$~~ & ~~$8.2\rightarrow8.3$ \\ \midrule
1     & $-7/2$ & 7.33 & 23.815 & 16.489 \\
2     & $-5/2$ & 7.21 & 23.671 & 16.467 \\
3     & $-3/2$ & 7.08 & 23.527 & 16.455 \\
4     & $-1/2$ & 6.94 & 23.381 & 16.450 \\
5     & $+1/2$ & 6.80 & 23.235 & 16.450 \\
6     & $+3/2$ & 6.64 & 23.088 & 16.455 \\
7     & $+5/2$ & 6.48 & 22.941 & 16.467 \\
8     & $+7/2$ & 6.31 & 22.794 & 16.489 \\ \bottomrule
\end{tabular}
\caption{HF split transition frequencies for CF level transitions  $8.1\rightarrow8.2$ ($x=0.1\%$), $8.1\rightarrow8.3$ ($x=0.01\%$) and $8.2\rightarrow8.3$ ($x=0.01\%$). All data are given in units of $\textrm{cm}^{-1}$. Uncertainties correspond to $\le\pm0.01$, $\le\pm0.001$ and $\le\pm0.003$\,$\textrm{cm}^{-1}$ \ for the $8.1\rightarrow8.2$, $8.1\rightarrow8.3$ and $8.2\rightarrow8.3$ transitions, respectively.}
\label{etab}
\end{table}

\section{Extraction of crystal field parameters and hyperfine interactions}
\label{sec: CF parameters}
The CF parameters of \lirefx{Ho} have been estimated previously based on CF level energies obtained as an average over their HF structure due to the limited resolution~\cite{gifeisman1978,christensen1979,goerllerwalrand1993,shakurov2005,babkevich2015}, or by magnetic susceptibility measurements~\cite{hansen1975,beauvillain1980,ronnow2007}. We improve on those earlier results by including the individually-resolved HF energies of all three CF transitions reported here and supplement these data with results from higher-lying CF states from Ref.~\cite{matmon2016}. We fit the CF parameters and the HF coupling constant $A_J$ simultaneously by numerically calculating the transition energies from the CF Hamiltonian (without HF interaction), as well as the HF splitting to first order in $A_J$. The transition energies are weighted with their measurement errors. This procedure only neglects small corrections to the CF energies due to HF interactions and the linear HF shift in $m_z$ due to second order terms in $A_J$, see Eqs.~(\ref{eq: HF corr doublet},\ref{eq: HF corr singlet}). The refined CF parameters are reported in Table~\ref{tab: cfp}, and we extract the HF coupling constant $A_J=0.02703\pm0.00003$\,$\textrm{cm}^{-1}$ \ in agreement with previous estimates in the literature of $A_J=0.0282\pm0.0005$~\cite{magarino1980} and $0.0270\pm0.0003$\,$\textrm{cm}^{-1}$ \ \cite{mennega1984}. The error bars of the CF parameters and $A_J$ are computed from the covariance matrix.

\begin{table}
\begin{tabular}{@{}l@{}|r@{}}
\toprule
CF parameter~ & Value [$\textrm{cm}^{-1}$]~~~~               \\ \midrule
~~~~~~$B_2^0$                  & $({-2.66}\pm0.05)\times10^{-1}$ \\
~~~~~~$B_4^0$                  & $({1.68}\pm0.04)\times10^{-3}$  \\
~~~~~~$B_4^{4}$                & $({2.81}\pm0.02)\times10^{-2}$ \\
~~~~~~$B_6^0$                  & $({5.74}\pm0.18)\times10^{-6}$  \\
~~~~~~$B_6^4$                  & $({5.60}\pm0.03)\times10^{-4}$  \\
~~~~~~$B_6^{-4}$               & $({0.00}\pm3.84)\times10^{-3}$  \\ \bottomrule
\end{tabular}
\caption{CF parameters extracted from the transition energy measurements. $B_4^{-4}$ is assumed to be zero~\cite{ronnow2007}.}
\label{tab: cfp}
\end{table}

A comparison of our CF parameter values in Table~\ref{tab: cfp} with previous results shows that we predict smaller values than previously, and we obtain significantly smaller values for $B_2^0$ and $B_4^0$. We attribute these corrections to the inclusion of the HF interaction term (to first order in $A_J$) in the Hamiltonian. In particular, fitting the HF structure allows us to use the magnetic moment of the 8.1 and 8.6 doublets (measured in Ref.~\cite{matmon2016}) as an additional constraint on the CF parameters, which determines the first order HF splitting. With the derived CF parameters we find a considerably ($\sim10 \%$) smaller magnetic moment $\mu/\mu_B= \pm 4.49$ of the 8.6 states than with previous CF parameters. We show the computed CF energies of the $^5$\textit{I}$_8$ manifold and their magnetic moments in Table~\ref{tab: elevels}. Compared to earlier reports, we find $2$-$10\%$ deviations for the predicted energies of the CF levels 8.7 to 8.13.

\begin{table}
\begin{tabular}{@{}l|c l c@{}}
\toprule
CF state~~ & ~Energy~[$\textrm{cm}^{-1}$]~& ~Symmetry  & ~~~$\braket{J_z} = \mu/(g_J \mu_\mathrm{B})$ \\ \midrule
~~~~~8.1  & 0~~~~    &  ~~~~~~$\Gamma_{3,4}$  &  ~\phantom{-}5.40       \\
~~~~~8.2  & 6.84   &  ~~~~~~$\Gamma_{2}$  &        \\
~~~~~8.3  & 23.31   &  ~~~~~~$\Gamma_{2}$  &        \\
~~~~~8.4  & 47.60  &  ~~~~~~$\Gamma_{1}$   &       \\
~~~~~8.5  & 56.92   &  ~~~~~~$\Gamma_{1}$   &       \\
~~~~~8.6  & 72.10 &  ~~~~~~$\Gamma_{3,4}$    &  ~-3.59    \\
~~~~~8.7  & 190.88~ &  ~~~~~~$\Gamma_{1}$   &       \\
~~~~~8.8  & 257.47~ &  ~~~~~~$\Gamma_{3,4}$  &   ~-2.30     \\
~~~~~8.9  & 275.31~ &  ~~~~~~$\Gamma_{2}$ &         \\
~~~~~8.10 & 275.38~ &  ~~~~~~$\Gamma_{1}$  &       \\
~~~~~8.11 & 288.66~ &  ~~~~~~$\Gamma_{1}$   &       \\
~~~~~8.12 & 294.65~  &  ~~~~~~$\Gamma_{3,4}$  &  ~\phantom{-}4.51      \\
~~~~~8.13 & 303.37~ &  ~~~~~~$\Gamma_{2}$   &      \\\bottomrule
\end{tabular}
\caption{Calculated energy levels of \ion{Ho} in \lirefx{Ho} \ based on the CF parameters shown in Table~\ref{tab: cfp}. The last column shows the expectation value $J_z$ of the magnetic $\Gamma_4$ states, which is proportional to their longitudinal magnetic moment $\mu$ and the Land\'e $g$-factor $g_J=5/4$.}
\label{tab: elevels}
\end{table}

\begin{figure}
\centering
\includegraphics[width=1\columnwidth]{Figures/2ndorder.pdf}
\caption{\noindent Energy differences $D_{2}$ (red) and $D_3$ (green) between neighboring HF transitions, shown as a function of the HF index. Linear fits to the experimental data are shown in orange. The respective slopes $s_{2,3}$ are a measure of the HF corrections $\propto m_z^2$. The numerical calculations based on our fitted CF parameters are shown in violet and cyan. The HF transition energies are shown in the inset for the $8.1\rightarrow8.2$ and $8.1\rightarrow8.3$ transition in the upper and lower panel, respectively.}
\label{2ndOrderFig}
\end{figure}

In contrast to $A_J$, the determination of the quadrupolar HF interaction constant $B$ requires precise knowledge of the deviations from the linear dipolar HF contributions. We utilize our high-resolution spectra ($8.1\rightarrow8.2$, $8.1\rightarrow8.3$, $8.2\rightarrow8.3$) to fit $B$ separately, by using the determined CF parameters and $A_J$, and numerically calculating the full HF spectrum. We find $B= 0.04 \pm 0.01$\,$\textrm{cm}^{-1}$, which is comparable to the literature value $B= 0.059$\,$\textrm{cm}^{-1}$ \ calculated from the free Ho atom~\cite{elliot1972}.

The parameters $A_J$ and $B$ allow us to numerically compute the HF spectrum. We present a comparison to the experimental data in Fig.~\ref{2ndOrderFig}. To emphasize the HF corrections $\propto m_z^2$, we look at the difference $D_n$ of transition frequencies $8.1\rightarrow8.n$ between neighboring $m_z$ for $n=2,3$. From Eqs.(\ref{eq: HF corr doublet},~\ref{eq: HF corr singlet}), we expect $D_{2,3}$ to be linear in $m_z$, with the slopes being a measure of the HF corrections $\propto m_z^2$. We find the slopes of $D_{2,3}$ (orange lines) to be $s_2=(7.2\pm0.5)\times10^{-3}$\,$\textrm{cm}^{-1}$ \ and $s_3=(6\pm1)\times10^{-4}$\,$\textrm{cm}^{-1}$, respectively, based on a linear regression. The numerically calculated values are shown in violet ($D_{2}$) and cyan ($D_{3}$). On average, we find the deviations of the experimental and numerically calculated values to be $16\%$ for the 8.2 and only $1.5\%$ for the 8.3 level, respectively. The error reflects the respective measurement resolutions.

The experimental values of $D_2$ and $D_3$ allow an order-of-magnitude estimation of the $m_z^2$-correction of the ground state (with the prefactor $\lambda_1/2$). We provide a detailed derivation thereof in the appendix~\ref{app: extraction of HF corr}. Namely, we neglect the quadrupolar interaction $B$ and restrict the sum over the CF states in Eqs.~(\ref{eq: HF corr doublet},~\ref{eq: HF corr singlet}) to the three lowest CF states, which contribute the most to the correction. We then exploit the anti-symmetry of the second-order corrections between the 8.1, 8.2 and 8.3 states to extract $\lambda_1=(s_1+s_2)/4= (2.0\pm0.2)\times10^{-3}$\,$\textrm{cm}^{-1}$ \ from our data, \cf \ Eq.~(\ref{eq: HF non-equidistance 1}). This is in agreement with the numerical calculation, yielding $\lambda_1=0.0024$\,$\textrm{cm}^{-1}$. Based on the errors found for the 8.2 and 8.3 energy level predictions, we expect a similar error of $\lesssim16\%$ for $\lambda_{1}$. Akin to $\lambda_1$, we estimate $\lambda_{2}=(-2.5\pm0.1)\times10^{-3}$\,$\textrm{cm}^{-1}$ \ and $\lambda_{3}=(1.9\pm0.3)\times10^{-3}$\,$\textrm{cm}^{-1}$ \ by also including the $8.2\rightarrow8.3$ transition. Both values are also in agreement with the numerical results $\lambda_{2}=-0.0040$\,$\textrm{cm}^{-1}$ \ and $\lambda_{3}=0.0017$\,$\textrm{cm}^{-1}$.

These $m_z^2$-corrections, \ie \ $\lambda_n$, have direct implications on the possibility to unambiguously address HF states, \eg \ in the context of quantum information processing. Specifically, $m_z\to m_z+1$ transitions within an electronic CF state $\ket{8.n}$ can only be driven if the HF line width is smaller than $\lambda_n$, which is the frequency difference of neighboring $m_z\to m_z+1$ transitions. At first sight, the line widths of our spectra seem not to satisfy this criterion. However, a quantitative evaluation of contributions to the line width is necessary to assess whether a regime (temperature, Ho-doping, etc.) for a specific CF state exists, where unambiguous addressing of HF states is possible~\cite{beckert2020c}. This problem can be circumvented by driving protocols involving another excited doublet CF level. Nuclear states can then be manipulated via a first $m_z$-conserving transition to an excited doublet with a subsequent transition to the $m_z$+1 state in the original CF level. For such manipulations involving a nuclear spin flip, the transition energies with different $m_z$ differ already in their first order hyperfine correction and do not rely on the much smaller $m_z^2$-corrections $\lambda$. Our spectra show that this condition is indeed fulfilled.

\section{Conclusions}
\label{conclusion}
We have extended the characterization of the ground CF state manifold of \lirefx{Ho} ($x=1\%$, $0.1\%$, and $0.01\%$) by direct optical measurements of transitions within the lowest three CF states. From the data we calculate the CF parameters, which differ from previous estimates because our refinement also considers the magnetic moments of the CF states as an additional fit constraint via the first order HF shift in $A_J$. In addition, this enables deducing the dipolar HF constant ${A_J=0.02703\pm0.00003}$\,$\textrm{cm}^{-1}$ \ purely by optical means. Using the CF parameters we predict the energies for the CF states of the $^5$\textit{I}$_8$ ground state manifold.
Our high measurement resolution allows us to determine the quadrupolar HF constant $B= 0.04\pm0.01$\,$\textrm{cm}^{-1}$ and subsequently to calculate the HF corrections of the three lowest CF states. We directly corroborate these calculations via estimations from our data. In addition, we report in appendix~\ref{Sec: Refractive Index} the far-infrared refractive index of \lirefx{Ho}. We conclude that specific addressing of individual HF transitions between doublet states is possible in \lirefx{Ho}, which is important in view of quantum information processing applications~\cite{grimm2020}.

\section{Acknowledgments}
FTIR spectroscopy data was taken at the X01DC beamline of the Swiss Light Source, Paul Scherrer Institut, Villigen, Switzerland. We thank H.~M.~R{\o}nnow, P. Babkevich and J.~Bailey for helpful discussions and experimental support. We thank S.~Stutz for technical support at the X01DC beamline. We acknowledge financial support by the Swiss National Science Foundation, Grant No. 200021\_166271, the European Research Council under the European Union's Horizon 2020 research and innovation programme HERO (Grant agreement No. 810451), and the Engineering and Physical Sciences Research Council, U.K. (
`HyperTerahertz' EP/P021859/1 and `COTS' EP/J017671/1).


text:
\section{Introduction}

Recent experiments on the electroproduction of $\pi^+$ at JLab
has drawn attention because the data from the experiment cover a wide
range of the photon momentum transfer $Q^2$ at high energy $W$
\cite{Blok:2008jy,Horn:2007ug,park}. Analyses of electroproduction
data have shown that the cross sections of the process are largely
determined by the pion and the nucleon electromagnetic form factors
that include information about the hadron structure
\cite{vgl,kaskulov,vrancx}. Therefore, electroproduction plays a
role not only in understanding the production mechanism  but also in
testing various sorts of form factors originating either from
phenomenological or from some theoretical basis on Quantum Chromodynamics (QCD)
\cite{kaskulov,gari1}.

One interesting feature of the process is that empirical data at
large momentum transfer $Q^2$ and high energy $W$ require a large
contribution of the $s$-channel proton pole to the transverse
cross section $d\sigma_T/dt$ in order to reproduce a set of four
separated differential cross sections,

    \begin{eqnarray}\label{crss1}
    &&2\pi\frac{d \sigma}{dt d\phi}= \frac{d\sigma_T}{dt}
    +\epsilon \frac{d\sigma_L}{dt} +\epsilon
    \frac{d\sigma_{TT}}{dt} \cos{2\phi}
    \nonumber\\&&

    +\sqrt{2\epsilon(\epsilon+1)}\frac{d\sigma_{LT}}{dt}
    \cos{\phi}\,,

    \end{eqnarray}

in a consistent manner. The longitudinal cross section
$d\sigma_L/dt$ follows the pion-pole dominance with the form
factor $F_\pi(Q^2)$ of a monopole type  from $\rho$-meson
dominance, as expected from
\begin{eqnarray}
d\sigma_L\propto \left|g_{\pi NN}F_\pi(Q^2)\right|^2.
\end{eqnarray}
For an enhancement of  $s$-channel exchange, therefore, in the
Regge model for $t$-channel meson exchange, Kaskulov and Mosel
(KM) \cite{kaskulov} introduced an $s$-dependence of the charge form
factor $F_s$ to the proton pole term to implement the contributions of
the $N^*$ resonances based on the Bloom-Gilman duality \cite{bloom},
i.e.,
\begin{eqnarray}\label{kmff}
F_s(Q^2,s)={\int_{M_p^2}^\infty ds_i{s_i^{-\beta}\over
s-s_i+i0^+}\left(1+\xi{Q^2\over s_i}\right)^{-2}\over
\int_{M_p^2}^\infty ds_i{s_i^{-\beta}\over s-s_i+i0^+}}\ ,
\end{eqnarray}
which results in Eq. (43) of Ref. 5.
Within a framework similar to that of KM, Vrancx and Ryckebusch (VR)
showed a fit of cross sections by using the proton charge form factor
$F_s$ of a dipole type  \cite{vrancx},
\begin{eqnarray}\label{vrff}
F_p(Q^2,s) = \left(1 + Q^2/\Lambda^2_{\gamma pp^{*}}(s)\right)^{-2},
\end{eqnarray}
in which case the cutoff mass has an energy $s$-dependence.
As shown in Ref.~6, however, the validity of such form
factors of proton in these works are questionable for the large $-t$
region, which is particularly true for  very recent data \cite{park}.
\\

In this work we investigate the process
\begin{eqnarray}\label{process}
\gamma^*(k)+ p(p)\to\pi^+(q)+ n(p')
\end{eqnarray}
based on an extension of our previous works \cite{bgyu1} to
electroproduction. Our purpose here is to find another set of
possible parameters for the pion and the proton charge form factors to fit
the cross sections. Avoiding complications in the parameterization
of the proton charge form factor as in Eqs. (\ref{kmff}) and
(\ref{vrff}), we reconsider the proton form factor of a simple
dipole type as in Eq. (\ref{vrff}), but now it has a constant cutoff
mass  that can be adjusted.  We then examine whether that proton form factor is valid up to
$-t\approx 5$ GeV$^2$ at large $Q^2$ and high energy $W$.

\begin{figure}[]
    \centering
    \includegraphics[width=0.9\columnwidth]{effdpv4.eps}
    \caption{(Color online) The $Q^2$-dependence of
    proton charge form factors used in the models of Refs. 5 and 6
    and the present work. The curves from KM and VR are estimated at
    $s=M_p^2$.
   }
    \label{emffs}
\end{figure}

For clarity, we work with a simple model consisting of $\pi+\rho$
Regge pole exchanges. In order to maintain gauge invariance of the
production amplitude, we have used the Gross-Riska prescription \cite{gross}
for proton and  pion charge form factors to constrain
the respective longitudinal components from coupling to virtual photon.
Thus, the production amplitude can be written as
\begin{eqnarray}\label{born}
\mathcal{M}=\bar{u}_N(p')\sqrt{2}\left[{\cal M}_{s,p} +{\cal
M}_{t,\pi} +{\cal M}_{t,\rho} \right] u_N(p),
\end{eqnarray}
with
\begin{eqnarray}\label{s}
    &&i{\mathcal M}_{s,p}=e g_{\pi NN}
    \gamma_5 {(\rlap{/}{p}+\rlap{/}{k}+M_p)\over s-M_p^2}
     \widetilde{F}_{1}(Q^2)\rlap{/}\epsilon\,,\\
    &&i{\mathcal M}_{t,\pi}=e g_{\pi NN}
    \widetilde{F}_{\pi}(Q^2){(2q-k)\cdot\epsilon\over t-m^2_\pi}\,\gamma_5 ,\\
    &&i{\mathcal M}_{t,\rho}=-i g_{\gamma\rho\pi} g_{\rho NN}
    F_{\rho}(Q^2)\epsilon^{\mu\nu\alpha\beta}{\epsilon_\mu k_\nu q'_\alpha \over
    t-m^2_\rho} \nonumber\\&&\hspace{1cm}\times
    \left( \gamma_\beta+\frac{\kappa_\rho}{4m_p}[\gamma_\beta,\rlap{/}{q'}]
    \right),
\end{eqnarray}
where $\widetilde{F}_{1}(Q^2)\rlap{/}{\epsilon}$ and
$\widetilde{F}_{\pi}(Q^2)(2q-k)^\mu$ are
\begin{eqnarray}
&&(F_1(k^2)-F_1(0))
\left(\rlap{/}\epsilon-\rlap{/}{k}{k\cdot\epsilon\over k^2}\right)
+ F_1(0)\rlap{/}\epsilon,\label{gr-1}
\\
&&(F_{\pi}(k^2)-F_{\pi}(0)) (2q-k)\cdot\left(\epsilon-
k{k\cdot\epsilon\over k^2}\right)
\nonumber\\&&\hspace{1cm}
+F_{\pi}(0)(2q-k)\cdot\epsilon,
\label{gr-2}
\end{eqnarray}
respectively. $Q^2=-k^2$ is the virtual photon momentum and $q'=q-k$
is the $t$-channel momentum transfer.  The gauge-invariant $\rho$
meson exchange is denoted by ${\cal M}_{t,\rho}$, with the
transition form factor $F_\rho(Q^2)$ at the $\gamma\rho\pi$
coupling vertex with the coupling constant
$g_{\gamma\rho\pi}=0.22$ GeV$^{-1}$. For the meson-baryon coupling
constants, we use $g_{\pi NN}=13.4$ and $g_{\rho NN}=3.4$ with
$\kappa_\rho=6.1$.
\\
\begin{figure}[]
    \centering
    \includegraphics[width=0.9\columnwidth]{dsdtdesy21dpv25.eps}
    \caption{(Color online) Differential cross sections at very forward
    angles $-t<0.3$ GeV$^2$ and small $Q^2$.
    The solid line results from the present work with $\Lambda_\pi=0.78$ GeV and
    $\Lambda_1=1.55$ GeV.
    The dotted line is from the KM model with $\Lambda_\pi=0.775$ (left) and 0.63 (right) GeV
    and $F_s(Q^{2},s)$ in Eq. (43)
    in Ref. 5.
    The dashed line is from the VR model with $\Lambda_\pi=0.655$ GeV and
    $F_p(Q^2,s)$ in Eq. (\ref{vrff}).
   Data are taken from Ref. 14
   for $Q^2=0.35$ GeV$^2$ and from Ref. 15
    for $Q^2=0.7$~GeV$^2$ (DESY).
} \label{desy21}
\end{figure}

To reggeize the $t$-channel meson exchange, we simply replace the
Feynmann pole with the Regge pole which can be collectively written
as
\begin{eqnarray}
&&{\cal R}^\varphi(s,t)\nonumber\\
&&=\frac{\pi\alpha'_\varphi}{\Gamma(\alpha_\varphi(t)+1-J)}
\frac{{1\over2}((-1)^J+e^{-i\pi\alpha_\varphi})}{\sin\pi\alpha_\varphi(t)}
\left(\frac{s}{s_0}\right)^{\alpha_\varphi(t)-J}
\end{eqnarray}
for the $\varphi$ meson of arbitrary spin $J$. In the present work,
we use the complex phase $e^{-i\pi\alpha(t)}$ for both the $\pi$ and
the $\rho$ exchanges with the trajectories
\begin{eqnarray}
\alpha_\pi (t) &=& 0.7(t-m_\pi^2)\,, \\
\alpha_\rho (t) &=& 0.83t+0.53\,,
\end{eqnarray}
respectively.
\\

\begin{figure*}[]
    \centering
    \includegraphics[width=1.33\columnwidth]{dsdt3expv16.eps}
    \caption{(Color online) Differential cross sections at forward angles $-t<1$ and
    $Q^2$ up to 4 GeV$^2$.
 The solid line results from the present work with $\Lambda_\pi=0.65$ GeV and $\Lambda=1.55$ GeV.
   Notations are the same as in Fig. \ref{desy21} for the dotted line with $\Lambda_\pi=0.68$ GeV
   and the dashed line with $\Lambda_\pi=0.655$ GeV.
   Data are taken from
   the F$\pi$-1 \cite{tadevosyan} (squares),
   F$\pi$-2 \cite{horn,Blok:2008jy} (triangles), and
   $\pi$-CT \cite{Horn:2007ug} (circles) experiments at JLab.}\label{jlab3exp}
\end{figure*}

Let us now discuss how to choose the pion and the proton charge form
factors in Eqs. (\ref{gr-1}) and (\ref{gr-2}), which are
needed to model the $\pi^+$ electroproduction process. In
the small $-t$ and low $Q^2$ region (see Fig. \ref{desy21} below,
for instance), cross section data show the longitudinal cross
section $d\sigma_L/dt$ to be large in comparison to the others. This is a
manifestation of the dominance of the pion exchange with the
charge form factor
\begin{eqnarray}\label{piff0}
F_{\pi}(Q^2)=\frac{1}{1+Q^2/\Lambda_\pi^2}\ ,
\end{eqnarray}
which is parameterized as a monopole-type from the vector meson
dominance. The cutoff mass is, therefore,  $\Lambda_\pi=m_\rho$,
which is somewhat larger than $\Lambda_\pi=0.71$ GeV fitted to the
measurement of the on-shell form factor in the elastic $e$N scattering
process. In the $\gamma^*p\to \pi^+n$ process, however, the pion
exchange proceeds via off-shell propagation; we, thus, consider
$\Lambda_\pi$ to be a fitting parameter to be varied in the range 0.6$\sim$
0.8 GeV  in this work.

\begin{figure*}[]
\centering
    \includegraphics[width=1.33\columnwidth]{dsudtjlabv11.eps}
    \caption{(Color online)
    Contribution of meson exchange in the unseparated cross section $d\sigma_U /dt$.
    Dash-dot-doted line is the contribution from $\pi$ exchange.
    Dash-dotted line corresponds to the contribution of
    $\pi+\rho$. Solid line results from $p({\rm proton})+\pi + \rho$ exchanges
    with $\Lambda_\pi=0.65$ GeV, and $\Lambda_1=1.55$ GeV.
    Data are taken from Ref. \cite{park}.} \label{clasdsdtu}
\end{figure*}

The proton exchange in the $s$-channel, though added to restore
gauge invariance, plays a role in the transverse cross section
$d\sigma_T/dt$, which corresponds to photoproduction of $\pi^+$ at
the photon point $Q^2=0$. It is reasonable to assume $N^*$
contributions in the resonance region \cite{kaskulov,vrancx}. The
$Q^2$ dependence of the proton form factor, called the Dirac form factor,
$F_1(Q^2)$
is determined from the Sachs electric and magnetic
form factors $G_E(Q^2)$ and $G_M(Q^2)$ by using the relation
\begin{eqnarray}\label{diracff}
F_1 (Q^2)=\frac{G_E(Q^2)+\tau G_M(Q^2)}{1+\tau},
\end{eqnarray}
where $\tau (Q^2)=Q^2/4M_P^2$ and $F_1(0)=1$. In the measurement
of  the on-shell form factors of a nucleon, $G_E(Q^2)$ and $G_M(Q^2)$
are applied and parameterized as
\begin{eqnarray}\label{sachff}
G_E(Q^2)=G_D(Q^2) \,, \,\,\,\, G_M(Q^2)=\mu_p G_D(Q^2),
\end{eqnarray}
where
\begin{eqnarray}\label{dpff}
G_D(Q^2)=\left(1+Q^2/\Lambda_1^2 \right)^{-2}
\end{eqnarray}
of the dipole type with $\Lambda_1^2 =0.71$ GeV$^2$
fitted to the empirical data \cite{hand}. $G_{E}(0)=1$ and
$G_{M}(0)=\mu_p=2.793$ normalized for the proton state. As
mentioned in the beginning, we adopt in this work the proton charge
form factor $F_1(Q^2)$ given in Eq. (\ref{diracff}), together with
Eqs. (\ref{sachff}), and (\ref{dpff}) above. Note that by
the definition in Eq. (\ref{diracff}), followed by Eqs. (\ref{sachff}) and
(\ref{dpff}), the proton charge form factor $F_1(Q^2)$ in this work
differs not only by the cutoff mass $\Lambda_{\gamma pp^*}(s)$
but also by the overall factor $\left(1+\tau\mu\over
1+\tau\right)$ from that of VR \cite{vrancx} given in Eq.
(\ref{vrff}).

In Fig. \ref{emffs}, we present the $Q^2$-dependence of the proton charge
form factors $F_1(Q^2)$ adopted in the models of KM, VR, and the
present work for comparison. The dashed line describes the form factor
$F_1(Q^2)$ in Eq. (\ref{dpff}) with $\Lambda_1 = 0.84$ GeV, which
is in good agreement with experimental data of Refs.
12 and 13.
The dotted line  is from the transition form
factor $F_s(Q^2, s\rightarrow M_p^2)$  in the KM model (Eq. (43) of
Ref 5)
while the dash-dotted line is from the form
factor in Eq. (\ref{vrff}) for the VR model. The solid line results
from the Dirac form factor $F_1(Q^2)$ with $\Lambda_1 = 1.55$ GeV used
in the present work.

While differences in the
proton form factor among the models are apparent, the cutoff mass
$\Lambda_\pi$ based on the pion form factor in Eq. (\ref{piff0}),
common to all models, should be different from one another for an
agreement with experiment. The KM model divides the cutoff mass into three
parts: $\Lambda_\pi=0.775$ GeV for $Q^2<0.4$ GeV$^2$, 0.63 GeV for
$0.6<Q^2<1.5$ GeV$^2$, and 0.68 GeV for $Q^2> 1.5$ GeV$^2$. The VR
model fixed $\Lambda_\pi=0.655$ GeV, regardless of the $Q^2$ range with a
smaller coupling constant $g_{\pi NN}=13.0$. We use
$\Lambda_\pi=0.78$ GeV for the DESY data and $\Lambda_\pi=0.65$ GeV
for the JLab data with the coupling constants concerning $\pi$
and $\rho$ being the same as those of the KM model.

Figure \ref{desy21} shows four differential cross sections
resulting from the KM model, the VR model, and the model used in this work for the small $-t$ and $Q^2$
region. Note that these model predictions are made within our
simple framework of ${\rm proton}+\pi+\rho$ exchanges where neglecting
the higher-spin meson exchange employed in the KM and the VR models, which is
by two orders of magnitude smaller than $\pi+\rho$. Therefore, the
differences between the models are basically the differences of the proton
form factors together with the cutoff $\Lambda_\pi$, as discussed
above. In this analysis, we fix the cutoff mass $\Lambda_1=1.55$
GeV over whole range of $Q^2$ and $-t$ for the proton form
factor in Eq. (\ref{diracff}) while the VM model uses $\Lambda_{\gamma
pp^*}(s)=0.84+(2.19-0.84)(1-M_p^2/s)$(GeV). The KM form factor
would yield $\Lambda_1\approx 1.48$ GeV in the limit of $s\to
M_p^2$, as estimated in correspondence with the same form factor
for the VM model \cite{vrancx}. Thus, all the cutoff masses for the proton form
factors in these models are almost two times larger than the
$\Lambda_1=0.84$ GeV fitted to the on-shell form factor from
experimental data, suggestive of  an $N^*$ contribution or an off-shell
effect in the $s$-channel.

\begin{figure}[]
\centering
      \includegraphics[width=0.83\columnwidth]{dsudtjlabmc.eps}
    \caption{(Color online)  $-t$ dependence of unseparated cross section $d\sigma_U /dt$.
Notations are the same as in Fig. \ref{jlab3exp} for  solid,
dotted and  dashed lines. The dash-dash-dotted line is from the VR
model with the inclusion of the strong hadronic form factor given
in Eq. (28) of Ref.\cite{vrancx}. Data are taken from Ref.
\cite{park}. }\label{clasdsdtumc}
\end{figure}

In the high-$Q^2$ region, the JLab data for small $-t$ and large $-t$
are presented in Fig. \ref{jlab3exp} and Fig. \ref{clasdsdtu}. Our
model reproduces those cross sections with a reduced value
$\Lambda_\pi=0.65$ GeV while using an fixed $\Lambda_p=1.55$ GeV, which is
independent of $Q^2$ and $-t$ momentum transfer.

The change in $\Lambda_\pi$ due to a variation of $Q^2$ can be
understood as a change in the measurement of the charge radius of pion,
\begin{eqnarray}
<r_\pi^2>=-6 \frac{dF_\pi}{dQ^2}\bigg|_{Q^2=0} =6/\Lambda_\pi^2,
\end{eqnarray}
which reads $<r_\pi^2>=0.384$ fm$^2$ for $\Lambda_\pi=0.78$ GeV
and $<r_\pi^2>=0.55$ fm$^2$ for $\Lambda_\pi=0.65$ GeV. This is
sensible because the size of the pion will increase with
higher resolution as $Q^2$ increases.

Based on the present model with cutoff masses consistent with DESY
and JLab data, we analyze the unseparated cross section
\begin{eqnarray}
d\sigma_U=d\sigma_T+\epsilon d\sigma_L
\end{eqnarray}
at the high $Q^2$ and large $-t$ region. We note in Fig.
\ref{clasdsdtu} that the contribution of $\rho$ exchange becomes
noticeable in this kinematical region. We should point out
that the sign of $\gamma\rho\pi$ coupling is of importance in
reproducing the data. A comparison of model predictions is given  in
Fig. \ref{clasdsdtumc} for the cross section at high $Q^2$ and
$-t$ up to 5 GeV$^2$. The result confirms the present work to be
more reliable for wider applications than the others.

\section*{Acknowledgment}
This work was partly supported  by grant NRF-2013R1A1A2010504
from the National Research Foundation (NRF) of Korea.