# Excitation of forbidden electronic transitions in atoms by Hermite-Gaussian
modes

Anton A. Peshkov111Corresponding author E-mail<EMAIL_ADDRESS>11 2
2 Elena Jordan 11 Markus Kromrey 11 Karan K. Mehta 33 Tanja E. Mehlstäubler 11
4 4 5 5 Andrey Surzhykov 11 2 2 6 6 Physikalisch-Technische Bundesanstalt,
D-38116 Braunschweig, Germany  Institut für Mathematische Physik, Technische
Universität Braunschweig, D-38106 Braunschweig, Germany  School of Electrical
and Computer Engineering, Cornell University, Ithaca, NY 14850, USA  Institut
für Quantenoptik, Leibniz Universität Hannover, D-30167 Hanover, Germany
Laboratorium für Nano- und Quantenengineering, Leibniz Universität Hannover,
Schneiderberg 39, 30167 Hannover, Germany  Laboratory for Emerging
Nanometrology Braunschweig, D-38106 Braunschweig, Germany

###### Abstract

Photoexcitation of trapped ions by Hermite-Gaussian (HG) modes from guided
beam structures is proposed and investigated theoretically. In particular,
simple analytical expressions for the Rabi frequencies of induced atomic
transitions are derived that depend both on the parameters of HG beams and on
the geometry of an experiment. By using these general expressions, we
investigate the ${}^{2}S_{1/2}\to\;^{2}F_{7/2}$ electric octupole (E3)
transition in an Yb+ ion, localized in the low–intensity center of the HG10
and HG01 beams. We show how the corresponding Rabi frequency can be enhanced
by properly choosing the polarization of incident light and the orientation of
an external magnetic field, which defines the quantization axis of a target
ion. The calculations, performed for experimentally feasible beam parameters,
indicate that the achieved Rabi frequencies can be comparable or even higher
than those observed for the conventional Laguerre-Gaussian (LG) modes. Since
HG-like modes can be relatively straightforwardly generated with high purity
and stability from integrated photonics, our results suggest that they may
form a novel tool for investigating highly-forbidden atomic transitions.

##  1 Introduction

In the past years there has been increasing interest in the use of structured
light fields in microscopy [1], optical tweezers [2], manipulating cold
quantum gases [3], as well as classical and quantum communication systems [4,
5]. This interest is spurred by the unique properties of structured light such
as spatially dependent amplitude, phase and polarization [6, 7]. From a
metrological perspective, structured light is of particular importance because
of its potential to excite narrow-line clock transitions with simultaneous
suppression of the undesirable light shift. This can be achieved since
structured light beams exhibit highly inhomogeneous intensity profiles with a
dark (low–intensity) center. An atom or ion, placed in this center, is exposed
to an electric field of low strength but high gradient, which favours
non–dipole transitions but weakly perturbs atomic levels. Such an “excitation
in the darkness” has been successfully demonstrated in experiments with a
single trapped 40Ca+ [8] and an 171Yb+ ion [9]. In these experiments, ions
interacted with conventional Laguerre-Gaussian (LG) modes and the excitation
probability exhibited a strong dependence on both polarization and applied
magnetic field orientation [10, 11]. Similar studies have been performed with
optical standing waves [12, 13].

Even though both experiments with 40Ca+ and 171Yb+ ions have provided
important steps towards operating the dipole-forbidden transitions with
structured light, they suffered from imprecise localization of a target ion
with respect to the low-intensity laser beam center. This uncertainty in
target localization arises not only from the thermal motion of an atom in a
trap, but also from the lack of laser beam pointing stability. The latter
problem can be partially solved by using the integrated schemes in which the
light delivery optics is fabricated directly into the atom traps [14, 15, 16].
Such photonic integrated circuits are produced using lithographic techniques
that result in rectangular waveguide structures. The geometry of these
structures allows fairly straightforward excitation of free-space beams
approximating Hermite-Gaussian (HG) modes, which in contrast to Laguerre-
Gaussian ones do not possess the rotational symmetry. In this contribution,
therefore, we propose to use HG modes in photonic integrated setups for
precision spectroscopy on strongly forbidden (atomic clock) transitions. To
assess the feasibility of this approach, we investigate here how a single
trapped Yb+ ion interacts with Hermite-Gaussian light.

The paper is organized as follows. In Sec. 2 we briefly discuss the derivation
of amplitudes for the excitation of a single trapped atom by polarized HG
beams. Based on the obtained expressions we explore in Sec. 3 the
laser–induced ${}^{2}S_{1/2}(F=0)\to\;^{2}F_{7/2}(F=3,M=0)$ electric octupole
(E3) transition in 171Yb+ ion localized in the dark center of the beam.
Special attention is paid to the question of how the Rabi frequency
$\Omega^{\rm(HG)}$ of this transition is affected by the polarization of
incident Hermite-Gaussian light and the orientation of an external magnetic
field. Moreover, we present the results of calculations of $\Omega^{\rm(HG)}$
for realistic experimental parameters which indicate that HG modes can serve
as a valuable tool for studying dipole–forbidden transitions in an integrated
optical setup. Finally, Sec. 4 provides a conclusion and an outlook.

##  2 Theory

###  2.1 Bessel modes

Before we analyze the process of photoexcitation by HG modes, we first recall
the simpler problem of an atom interacting with Bessel light. Since this topic
has been extensively discussed in the literature [11, 19, 17, 18], we restrict
ourselves here to basic expressions. The Bessel beam is characterized by a
well-defined projection $m_{\gamma}$ of the total angular momentum upon its
propagation direction, the helicity $\lambda$, and the longitudinal component
$k_{z}$ of the linear momentum. Moreover, the absolute value of the transverse
momentum $|\bm{k}_{\perp}|=\varkappa$, and hence the frequency
$\omega=ck=c\sqrt{k_{z}^{2}+\varkappa^{2}}$, are also fixed. The vector
potential of such a Bessel beam

$\displaystyle\bm{A}^{(\text{B})}_{m_{\gamma}\lambda}(\bm{r};\,\varkappa,k_{z})=\int
a_{\varkappa
m_{\gamma}}(\bm{k}_{\perp})\,\bm{e}_{\bm{k}\lambda}\,e^{i\bm{k}\bm{r}}\,\frac{d^{2}\bm{k}_{\perp}}{(2\pi)^{2}}\,$
(1)

can be written as a superposition of plane waves with wave vectors $\bm{k}$
uniformly distributed upon the surface of a cone, whose axis coincides with
the light propagation direction and whose opening angle is given by
$\theta_{k}=\arcsin(\varkappa/k)$. In Eq. (1), $\bm{e}_{\bm{k}\lambda}$ is the
photon polarization vector satisfying the Coulomb gauge condition
$\bm{e}_{\bm{k}\lambda}\cdot\bm{k}=0$, and the amplitude $a_{\varkappa
m_{\gamma}}(\bm{k}_{\perp})$ is defined by:

$\displaystyle a_{\varkappa
m_{\gamma}}(\bm{k}_{\perp})=\frac{2\pi}{\varkappa}\,(-i)^{m_{\gamma}}\,e^{im_{\gamma}\phi_{k}}\,\delta(k_{\perp}-\varkappa)\,,$
(2)

where $\phi_{k}$ is the azimuthal angle of the wave vector.

With the help of the vector potential (1), we can write down the amplitude for
a radiative transition
$|\alpha_{g}F_{g}\,M_{g}\rangle+\gamma\to|\alpha_{e}F_{e}\,M_{e}\rangle$. The
ground and excited states of an atom are characterized here by nuclear $I$,
electron $J$, and total $\bm{F}=\bm{I}+\bm{J}$ angular momenta, and the
projection $M=M_{F}$ on the atomic quantization axis. This axis is directed
along the external magnetic field $\bm{B}$, tilted at an angle $\theta$ with
respect to the light propagation direction (see Fig. 1). Moreover, $\alpha$
denotes all the additional quantum numbers that are needed for a unique
specification of the states. The transition amplitude is:

$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{B})}$
$\displaystyle=\Bigg{\langle}\alpha_{e}F_{e}M_{e}\Bigg{|}\sum_{q}\bm{\alpha}_{q}\,\bm{A}_{m_{\gamma}\lambda}^{(\text{B})}(\bm{r}_{q};\,\varkappa,k_{z})\Bigg{|}\alpha_{g}F_{g}M_{g}\Bigg{\rangle}\,$
$\displaystyle=\int a_{\varkappa
m_{\gamma}}(\bm{k}_{\perp})\,e^{-i\bm{k}_{\perp}\bm{b}}\,$
$\displaystyle\times\Bigg{\langle}\alpha_{e}F_{e}M_{e}\Bigg{|}\sum_{q}\bm{\alpha}_{q}\,\bm{e}_{\bm{k}\lambda}\,e^{i\bm{k}\bm{r}_{q}}\Bigg{|}\alpha_{g}F_{g}M_{g}\Bigg{\rangle}\,\frac{d^{2}\bm{k}_{\perp}}{(2\pi)^{2}}\,,$
(3)

where $\bm{\alpha}_{q}$ is the vector of Dirac matrices for the $q$th
electron, and the impact parameter $\bm{b}$ specifies the lateral position of
an atom with respect to the beam axis. This parameter is introduced since
Bessel and Hermite-Gaussian modes have a complex position-dependent structure.
For Bessel beams, the case of zero impact parameter, $b=0$, corresponds to an
atom located on a vortex line [11].

Figure 1: Geometry of the excitation of a single trapped atom by
Hermite–Gaussian modes. The angle $\theta$ defines the direction of the
quantization axis (determined by the external magnetic field) with respect to
the light propagation direction. The $x$-$z$ plane is spanned by the light
propagation axis and the magnetic field. The atom is localized in the vicinity
of the beam center located at $x=y=0$.

In order to compute the amplitude (2.1), we need to choose the quantization
axis of the entire system “atom plus light”. The two most obvious choices are
the light propagation axis and the magnetic field direction. Of course,
physical observables are independent of a particular choice of coordinate
system. For the purposes of analysis, however, it is practical to take the
overall quantization axis along the $\bm{B}$-field. Having defined the
geometry, we can now apply the standard multipole decomposition technique to
the vector potential of the radiation field [20]. Namely, the expansion of the
plane-wave components of $\bm{A}^{(\text{B})}_{m_{\gamma}\lambda}$ is:

$\displaystyle\bm{e}_{\bm{k}\lambda}\,e^{i\bm{k}\bm{r}}=$
$\displaystyle\sqrt{2\pi}\sum_{pLMM^{\prime}}i^{L}\,[L]^{1/2}\,(i\lambda)^{p}\,D^{L}_{M\lambda}(\phi_{k},\,\theta_{k},0)\,$
$\displaystyle\times
D^{L}_{M^{\prime}M}(\pi,\,\theta,\,\pi)\,\bm{a}^{(p)}_{LM^{\prime}}(\bm{r})\,,$
(4)

where $D^{L}_{M\lambda}$ is the Wigner $D$-function, $[L]=2L+1$, and
$\bm{a}^{(p)}_{LM^{\prime}}$ refers to magnetic $(p=0)$ and electric $(p=1)$
multipole potentials. In Eq. (2.1), $\theta_{k}$ and $\phi_{k}$ are angles
that determine the propagation direction of each plane-wave component with
respect to the beam axis and the $x$-$z$ plane, respectively. Making use of
the Wigner-Eckart theorem and integrating over $\bm{k}_{\perp}$, we find:

$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{B})}(\bm{b})=\sqrt{2\pi}\sum_{pLM}i^{L}\,[L,F_{g}]^{1/2}\,(i\lambda)^{p}\,(-1)^{m_{\gamma}}\,i^{M}\,$
$\displaystyle\;\;\;\times
e^{i(m_{\gamma}-M)\phi_{b}}\,d^{L}_{M,\lambda}(\theta_{k})\,d^{L}_{M,M_{e}-M_{g}}(\theta)\,J_{m_{\gamma}-M}(\varkappa
b)\,$ $\displaystyle\;\;\;\times\langle
F_{g}M_{g}\,LM_{e}-M_{g}|F_{e}M_{e}\rangle\,(-1)^{J_{e}+I+F_{g}+L}\,$
$\displaystyle\;\;\;\times\left\\{\begin{array}[]{ccc}F_{e}&F_{g}&L\\\
J_{g}&J_{e}&I\end{array}\right\\}\,\langle\alpha_{e}J_{e}||H_{\gamma}(pL)||\alpha_{g}J_{g}\rangle\,,$
(7)

with $d^{L}_{M,\lambda}$ and $J_{m_{\gamma}-M}$ being the small Wigner and
Bessel functions, respectively. Here we have introduced, moreover, the reduced
matrix element of the transition operator
$H_{\gamma}(pL)=\sum_{q}\bm{\alpha}_{q}\bm{a}^{(p)}_{L}(\bm{r}_{q})$ that does
not depend on the projections of the angular momenta of an atom and photon. As
seen from Eq. (2.1), the matrix element
$\mathcal{M}_{M_{e}M_{g}}^{(\text{B})}$ depends on the tilt angle $\theta$ of
the magnetic field, the atom’s impact parameter $\bm{b}$, the opening angle
$\theta_{k}$, the helicity $\lambda$, and the total angular momentum
projection $m_{\gamma}$ of the beam. In the past, Eq. (2.1) has been employed
to study transitions in atoms and ions exposed to Bessel and Laguerre-Gaussian
beams [9, 21]. The application of this formula to the latter case is justified
when an atom is located near a vortex line, where the phase and intensity
patterns of LG and Bessel modes are almost identical [22].

Figure 2: Top: Intensity profiles for Hermite–Gaussian modes HG00 (left), HG10
(middle), and HG01 (right), linearly polarized along the $x$-axis. Bottom:
Corresponding probabilities of the
$|F_{g}\\!=\\!0,\,M_{g}\\!=\\!0\rangle+\gamma\to\,|F_{e}\\!=\\!3,\,M_{e}\\!=\\!0\rangle$
transition as a function of the magnetic field tilt angle $\theta$ for a
single 171Yb+ ion placed exactly at the beam center (blue solid lines) or
exhibiting spatial distribution with width parameter $\sigma=40$ nm (red
dashed lines) and $\sigma=100$ nm (green dash-dotted lines). Calculations were
performed for the wavelength of 467 nm, the beam waist $w_{0}=3$ $\mu$m, the
total power $P=1$ mW, and the pulse duration $t=1.5$ ms.

###  2.2 Hermite-Gaussian modes

Having discussed how to compute the transition amplitudes for Bessel (and LG)
radiation, we can consider now photoexcitation of an atom by Hermite-Gaussian
modes HGmn. These are known to be solutions of the paraxial wave equation in
Cartesian coordinates and are usually characterized by two indices $m$ and
$n$, which determine their transverse intensity profile [23]. In our study we
will focus on the modes HG00 (Gaussian beam), HG10, and HG01 that are planned
to be used in future clock on a chip experiments at PTB.

In contrast to diffraction-free Bessel solutions, the intensity distribution
of which is the same in every plane normal to the beam axis, the width of HG
intensity profile changes during propagation. In planned experiments, a target
atom or ion will be placed in the plane of minimum width, at $z=0$. In this
plane the vector potentials of the HG00, HG10, and HG01 modes with the
frequency $\omega$ and the beam waist $w_{0}$ can be written as:

$\displaystyle\bm{A}^{(\text{HG}_{00})}(x,y)$
$\displaystyle=\bm{e}A_{0}\,e^{-(x^{2}+y^{2})/w_{0}^{2}}\,,$ (8a)
$\displaystyle\bm{A}^{(\text{HG}_{10})}(x,y)$
$\displaystyle=\bm{e}A_{0}\,\frac{2x}{w_{0}}\,e^{-(x^{2}+y^{2})/w_{0}^{2}}\,,$
(8b) $\displaystyle\bm{A}^{(\text{HG}_{01})}(x,y)$
$\displaystyle=\bm{e}A_{0}\,\frac{2y}{w_{0}}\,e^{-(x^{2}+y^{2})/w_{0}^{2}}\,,$
(8c)

where $\bm{e}$ is the polarization vector [24]. In Eqs. (8), the constant
$A_{0}=\sqrt{4P/(c\epsilon_{0}\omega^{2}\pi w_{0}^{2})}$ is chosen in such a
way that an integral of the intensity
$I=c\epsilon_{0}\omega^{2}|\bm{A}|^{2}/2$ over the beam cross-section gives
the total power $P$ [25, 26]. It follows from Eqs. (8) that the HG10 and HG01
modes have one nodal line directed along $y$ and $x$ axes, respectively. As
seen from the upper panel of Fig. 2, each of these lines separates two bright
spots. Thus, in contrast to the case of Bessel and Gaussian modes, the HG10
and HG01 solutions do not have axial symmetry with respect to the direction of
light propagation.

Figure 3: Same as Fig. 2, but for the linear polarization along the $y$-axis.

Rigorous calculation of the transition amplitudes for Hermite-Gaussian modes
is somewhat more complicated than for the case of Bessel radiation.
Fortunately, in the region near the beam center located at $x=y=0$ the HG00,
HG10, and HG01 modes (8) can be well approximated by a linear combination of
Bessel solutions (1) in the paraxial regime. In particular, for the circularly
polarized Hermite-Gaussian beams, we find:

$\displaystyle\bm{A}^{(\text{HG}_{00})}_{\lambda=\pm 1}$
$\displaystyle\approx\pm iA_{0}\,\bm{A}^{(\text{B,
par})}_{m_{l}=0,\,\lambda=\pm 1}\,,$ (9a)
$\displaystyle\bm{A}^{(\text{HG}_{10})}_{\lambda=\pm 1}$
$\displaystyle\approx\pm 1.04\,\frac{iA_{0}}{\sqrt{2}}\left[\bm{A}^{(\text{B,
par})}_{m_{l}=+1,\,\lambda=\pm 1}-\bm{A}^{(\text{B,
par})}_{m_{l}=-1,\,\lambda=\pm 1}\right]\,,$ (9b)
$\displaystyle\bm{A}^{(\text{HG}_{01})}_{\lambda=\pm 1}$
$\displaystyle\approx\pm 1.04\,\frac{A_{0}}{\sqrt{2}}\left[\bm{A}^{(\text{B,
par})}_{m_{l}=+1,\,\lambda=\pm 1}+\bm{A}^{(\text{B,
par})}_{m_{l}=-1,\,\lambda=\pm 1}\right]\,,$ (9c)

where

$\displaystyle\bm{A}_{m_{l}\lambda}^{(\text{B, par})}(\bm{r})$
$\displaystyle\approx\bm{A}_{m_{\gamma}\lambda}^{(\text{B})}(\bm{r};\,\varkappa\ll
k_{z})\,$ $\displaystyle=\bm{e}_{\lambda}\,(-i)^{\lambda}J_{m_{l}}(\varkappa
r_{\perp})\,e^{im_{l}\phi}\,e^{ik_{z}z}\,.$ (10)

The latter expression is obtained from Eq. (1) by integrating over
$\bm{k}_{\perp}$ under the paraxial condition $\varkappa\ll k_{z}$, see Ref.
[11] for further details. In Eq. (2.2), $m_{l}=m_{\gamma}-\lambda$ denotes the
projection of the orbital angular momentum, and
$\bm{e}_{\lambda}\equiv\bm{e}_{\bm{k}\parallel z,\lambda}$ with $\lambda=\pm
1$ stands for the vector describing the states of right-hand and left-hand
circularly polarized radiation.

In order to approximate the HG solutions (8) by the paraxial Bessel beams
(2.2), we have also adjusted the ratio of the transverse to the longitudinal
photon momentum components determined by the opening angle $\theta_{k}$, as
well as the general prefactors. In particular, while for HG00 we have chosen
$\theta_{k}=\arcsin[2/(w_{0}k)]$, the opening angle
$\theta_{k}=\arcsin[2.6/(w_{0}k)]$ and the prefactor 1.04 are taken for HG10
and HG01. For these choice of parameters, Eqs. (8) and (9) give almost
equivalent vector potentials in the vicinity of the beam center located at
$x=y=0$. For example, for the case of the modes with wavelength 467 nm and
beam waist 3 $\mu$m displayed in the upper panel of Fig. 2, this region is
bounded by $-2$ $\mu$m $<$ $x$ $<$ $2$ $\mu$m and $-2$ $\mu$m $<$ $y$ $<$ $2$
$\mu$m. This is sufficient to describe future experiments with Doppler cooled
atoms in which atom localization near the beam center below 40 nm is expected
due to thermal fluctuations. For the chosen parameters, we can also neglect
the effects of the $z$-dependence of the beam intensity profile, since its
Rayleigh range of about 60 $\mu$m (i.e. the distance at which the beam width
increases noticeably) is much larger than the expected fluctuations of the $z$
coordinate of an atom.

By making use of Eqs. (9), we can express the transition amplitudes for
excitation of an atom by circularly polarized HG modes

$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00})}(\lambda=\pm 1)$
$\displaystyle\approx\pm i\,A_{0}\,\mathcal{M}_{M_{e}M_{g}}^{(\text{B,
par})}(m_{l}=0,\,\lambda=\pm 1)\,,$ (11a)
$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{10})}(\lambda=\pm 1)$
$\displaystyle\approx\pm 1.04A_{0}\,\frac{i}{\sqrt{2}}$
$\displaystyle\times\left[\mathcal{M}_{M_{e}M_{g}}^{(\text{B,
par})}(m_{l}=+1,\,\lambda=\pm 1)\right.$
$\displaystyle\left.\;\;-\mathcal{M}_{M_{e}M_{g}}^{(\text{B,
par})}(m_{l}=-1,\,\lambda=\pm 1)\right]\,,$ (11b)
$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{01})}(\lambda=\pm 1)$
$\displaystyle\approx\pm 1.04A_{0}\,\frac{1}{\sqrt{2}}$
$\displaystyle\times\left[\mathcal{M}_{M_{e}M_{g}}^{(\text{B,
par})}(m_{l}=+1,\,\lambda=\pm 1)\right.$
$\displaystyle\left.\;\;+\mathcal{M}_{M_{e}M_{g}}^{(\text{B,
par})}(m_{l}=-1,\,\lambda=\pm 1)\right]\,,$ (11c)

in terms of their Bessel counterparts (2.1). These formulas can be readily
used to obtain the amplitudes for linearly polarized light:

$\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(x)=$
$\displaystyle\frac{1}{\sqrt{2}}\left[\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(\lambda=+1)\right.$
$\displaystyle\left.\;\;\;\;+\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(\lambda=-1)\right]\,,$
(12a) $\displaystyle\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(y)=$
$\displaystyle\frac{i}{\sqrt{2}}\left[\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(\lambda=-1)\right.$
$\displaystyle\left.\;\;\;\;-\mathcal{M}_{M_{e}M_{g}}^{(\text{HG}_{00,10,01})}(\lambda=+1)\right]\,,$
(12b)

where we have employed the standard relationships
$\bm{e}_{x}=(\bm{e}_{\lambda=+1}+\bm{e}_{\lambda=-1})/\sqrt{2}$ and
$\bm{e}_{y}=i(\bm{e}_{\lambda=-1}-\bm{e}_{\lambda=+1})/\sqrt{2}$ between
linear and circular polarization unit vectors [20].

##  3 Results and discussion

While the theory presented above can be applied to describe excitation of an
arbitrary atom by low-order HG modes, here we will consider the
$4f^{14}6s\;^{2}S_{1/2}(F\\!=\\!0)\to\;4f^{13}6s^{2}\;{}^{2}F_{7/2}(F\\!=\\!3)$
transition in the 171Yb+ ion. This 467 nm transition, proceeding via the E3
channel, is of great interest since it serves as the reference for high-
precision optical clocks and enables searches for new physics [27, 28, 29].
Below we shall focus our attention on the transition between the magnetic
sublevels $M_{g}=M_{e}=0$ that can be spectroscopically resolved by the Zeeman
effect. As seen from Eq. (2.1), the evaluation of the transition amplitudes
requires a knowledge of the reduced matrix element $\langle
4f^{13}6s^{2}\,{}^{2}F_{7/2}||H_{\gamma}(E3)||4f^{14}6s\,^{2}S_{1/2}\rangle$.
Its calculation for 171Yb+ is a very complicated task, requiring the use of
sophisticated atomic structure theories. To avoid this, we can estimate the
matrix element from the expression

$\displaystyle\langle
4f^{13}6s^{2}\,{}^{2}F_{7/2}||H_{\gamma}(E3)||4f^{14}6s\,^{2}S_{1/2}\rangle=\sqrt{\frac{1}{\pi\alpha\omega\tau}}\,,$
(13)

which relates it to the measured lifetime $\tau=4.98\times 10^{7}$ s of the
${}^{2}F_{7/2}$ excited state [30, 31]. We shall also assume that the total
power of the beam is $P=1$ mW and its waist is $w_{0}=3$ $\mu$m. A similar
parameter range has been employed in recent experiments on excitation of the
E3 transition in 171Yb+ ion by LG light [9].

###  3.1 Transition probability

We are ready now to use Eqs. (12) to compute the probability $W$ of excitation
of an 171Yb+ ion by linearly polarized HG modes. For a well-defined impact
parameter of an atom, this transition probability is given by
$W(\bm{b})=\left|ec\mathcal{M}^{(\text{HG})}(\bm{b})/\hbar\right|^{2}t^{2}/4$.
This approximate formula was derived from the analysis of Rabi oscillations
for short interaction times $t$ [32]. In calculations below we assume $t=1.5$
ms which is much smaller than the inverse Rabi frequency.

As already mentioned above, in real experiments it is not possible to achieve
perfect control of the atom’s position. To account for such uncertainty, we
assume that the impact parameter $\bm{b}$ follows a Gaussian distribution with
a width $\sigma$:

$\displaystyle
f(\bm{b})=\frac{1}{2\pi\sigma^{2}}e^{-\frac{\bm{b}^{2}}{2\sigma^{2}}}\,.$ (14)

With the help of $f(\bm{b})$, we can use a semi-classical approximation to
express the transition probability for a single atom target, centered on the
beam axis, as:

$\displaystyle W^{(\text{HG})}_{x,\,y}=\frac{(ect)^{2}}{4\hbar^{2}}\int
f(\bm{b})\,\left|\mathcal{M}_{M_{e}M_{g}}^{(\text{HG})}\left(x,y\right)\right|^{2}\,d^{2}\bm{b}\,.$
(15)

In the past this approximation has been successfully employed to describe
experimental results obtained for LG modes [9].

Fig. 2 shows the transition probabilities as a function of the magnetic field
tilt angle $\theta$ for HG00 (left), HG10 (middle), and HG01 (right) modes,
linearly polarized along the $x$-direction. Here, the blue solid line
represents the results for an atom placed exactly at the beam center ($b=0$,
$\sigma=0$), while the red dashed and green dash-dotted lines correspond to
the predictions of Eq. (15) for the spatial distribution $f(\bm{b})$ with
width parameters $\sigma=40$ nm and $100$ nm, respectively. Since beam-
pointing fluctuations are expected to be negligible in an integrated optics
setup, we estimated these parameters based on the thermal spread of an spatial
atomic wavepacket. In this case, $\sigma$ can be derived from the relation
$m_{\text{a}}\omega_{\text{r}}^{2}\sigma^{2}/2=k_{\text{B}}T/2$ for a
classical harmonic oscillator, where $m_{\text{a}}$ is the mass of an ion and
$T$ is its temperature. With the help of this expression, we obtained the
above mentioned widths $\sigma$ for a realistic trapping frequency
$\omega_{\text{r}}=2\pi\times 600$ kHz and temperatures
$T=T_{\text{Doppler}}=0.5$ mK corresponding to the Doppler cooling limit and
$T=3.1$ mK. In the later discussion of the light shift we will also consider a
smaller width $\sigma=10$ nm, which corresponds to the size of the quantum
mechanical ground-state wave function of the trapped Yb+ ion. We see that the
excitation probability $W^{(\text{HG}_{00})}_{x}$ for the Gaussian beam is
insensitive to variations of the target size. This can be expected since the
Gaussian beam size $\approx 3$ $\mu$m is large compared to any of the width
parameters considered here. Therefore, the intensity and phase of HG00 mode
are almost constant over the entire spread of a target. For this case, the
excitation probability exhibits the well-known plane-wave behavior [11],
reaching maximum values at angles $\theta=31^{\circ}$, $90^{\circ}$, and
$149^{\circ}$.

A different behavior of $W^{(\text{HG})}_{x}$ can be observed for HG10 and
HG01 modes. For these two cases, the delocalization of an atom may influence
the excitation probability; however, the effect becomes pronounced only for
very large targets, $\sigma\approx 100$ nm. Such a $\sigma$-dependence is
caused by the inhomogeneity of the electric field distribution of HG10 and
HG01 beams. Moreover, their complex internal structure also affects the
$\theta$-dependence of the transition probability. For example,
$W^{(\text{HG}_{10})}_{x}$ is enhanced at $\theta=0^{\circ}$, $60^{\circ}$,
$120^{\circ}$, and $180^{\circ}$, which is much different from what was
observed for the Gaussian regime. In contrast, $W^{(\text{HG}_{01})}_{x}$
exhibits a qualitatively similar $\theta$-dependence as
$W^{(\text{HG}_{00})}_{x}$, but is strongly suppressed and becomes observable
only for large targets.

Figure 4: Rabi frequencies (17) as a function of the magnetic field tilt angle
$\theta$ for HG00 (black line) and HG10 (blue line) modes linearly polarized
in the $x$-direction and for HG01 (red line) mode linearly polarized in the
$y$-direction. The atom is placed in the beam center ($b=0$, $\sigma=0$). All
other parameters are the same as in Fig. 2.

Fig. 3 displays the excitation probabilities for the same three modes, HG00,
HG10, and HG01, but for linear polarization in the $y$-direction. As seen from
the figure, HG00 and HG10 light beams hardly induce the
$M_{g}\\!=\\!0\to\,M_{e}\\!=\\!0$ transition, regardless of the tilt angle
$\theta$ of the magnetic field. This suppression of the E3 transition can be
explained using symmetry arguments. For example, for plane-wave radiation and
the transition of interest, a detailed discussion of zero excitation
probability based on the symmetry analysis of ionic states and of incident
light has been presented in Ref. [11]. Such an analysis can be extended to
HG00 and HG10 modes when an atom is placed on the vortex line, $b=0$, but is
not fully justified for a delocalized atom. Nevertheless, for the latter case,
$\sigma>0$, the probability of the E3 transition is very tiny and hence is not
seen in the left and middle panels of Fig. 3. In contrast, excitation with
HG01 is possible for small and large tilt angles and is most pronounced for
$\theta=0^{\circ}$ and $180^{\circ}$.

###  3.2 Rabi frequency and light shift

As seen from Figs. 2 and 3, the results obtained for a perfectly localized
atom, $b=0$ and $\sigma=0$, differ only little from the predictions of Eq.
(15) for $\sigma=40$ nm at Doppler temperature. We argue therefore that the
atomic spread effects can be neglected in the analysis of future experiments
in which the size of a target is much smaller than the beam waist. As
mentioned above, this is the case for planned studies at PTB, where an
integrated optics scheme improves the pointing stability of the laser beam
relative to the ion. The approximation of a perfectly localized atom allows
one to derive simple analytical expressions for the Rabi frequency:

$\displaystyle\Omega^{(\text{HG})}=\frac{ec}{\hbar}\left|\mathcal{M}_{M_{e}M_{g}}^{(\text{HG})}(b=0)\right|\,.$
(16)

Indeed, by making use of Eqs. (12) and (16) for the case of the
$|F_{g}\\!=\\!0,\,M_{g}\\!=\\!0\rangle+\gamma\to\,|F_{e}\\!=\\!3,\,M_{e}\\!=\\!0\rangle$
$E3$ transition, induced by various linearly polarized Hermite-Gaussian modes,
we find:

$\displaystyle\Omega^{(\text{HG}_{00})}_{x}$
$\displaystyle\approx\sqrt{\frac{6P}{\hbar\omega
k^{2}w_{0}^{2}\tau}}\,\left|\sin\theta(5\cos^{2}\theta-1)\right|\,,$ (17a)
$\displaystyle\Omega^{(\text{HG}_{10})}_{x}$ $\displaystyle\approx
2.7\sqrt{\frac{3P}{\hbar\omega
k^{4}w_{0}^{4}\tau}}\,\left|\cos\theta(15\cos^{2}\theta-11)\right|\,,$ (17b)
$\displaystyle\Omega^{(\text{HG}_{01})}_{y}$ $\displaystyle\approx
2.7\sqrt{\frac{3P}{\hbar\omega
k^{4}w_{0}^{4}\tau}}\,\left|\cos\theta(5\cos^{2}\theta-1)\right|\,,$ (17c)
$\displaystyle\Omega^{(\text{HG}_{00})}_{y}$
$\displaystyle=\Omega^{(\text{HG}_{10})}_{y}=\Omega^{(\text{HG}_{01})}_{x}\approx
0\,,$ (17d)

where we have employed the paraxial approximation of the Bessel amplitude
obtained from Eq. (2.1) for small opening angles $\theta_{k}$.

In Fig. 4 we show non-zero Rabi frequencies $\Omega^{(\text{HG}_{00})}_{x}$,
$\Omega^{(\text{HG}_{10})}_{x}$, and $\Omega^{(\text{HG}_{01})}_{y}$,
calculated for the set of parameters mentioned at the beginning of Sec. 3. We
see that while the Rabi frequency may reach
$\Omega^{(\text{HG}_{00})}_{x}\approx 2\pi\times 90$ Hz for the Gaussian mode
and $\theta=31^{\circ}$ or $149^{\circ}$, its maximum values are
$\Omega^{(\text{HG}_{10})}_{x}=\Omega^{(\text{HG}_{01})}_{y}\approx 2\pi\times
12$ Hz for the HG10 and HG01 light beams and $\bm{B}\parallel\bm{k}_{z}$. It
is also informative to compare $\Omega^{(\text{HG}_{10})}$ and
$\Omega^{(\text{HG}_{01})}$ with the corresponding Rabi frequencies for
Laguerre-Gaussian modes that were used in the previous experiments [9]. By
using Eq. (2.2) and the theory from above, we find that
$\Omega^{(\text{HG}_{10})}_{x}/\Omega^{(\text{LG}_{01})}_{x}=\Omega^{(\text{HG}_{01})}_{y}/\Omega^{(\text{LG}_{01})}_{y}=\sqrt{2}$,
thus indicating that excitation probability of the E3 transition induced by HG
modes is two times larger than that obtained by LG light. Moreover, we
estimated for our beam parameters and an optimal excitation geometry that an
171Yb+ ion placed in the center of HG10 or HG01 modes experiences the light
shift in the range from $\Delta\nu=50$ Hz to 8 Hz if the target size varies
from $\sigma=100$ nm to 40 nm [30, 33]. An even smaller light shift of about
$\Delta\nu=0.5$ Hz can be achieved for the “quantum-mechanical” width
$\sigma=10$ nm. These values are comparable to light shifts for LG01 modes and
are much smaller than $\Delta\nu=205$ Hz obtained for a Gaussian beam with
power adjusted to provide a similar Rabi frequency.

##  4 Conclusions and outlook

In this paper we propose using Hermite-Gaussian modes in photonic integrated
circuits to induce dipole-forbidden atomic clock transitions. In order to
investigate the feasibility of this scheme, we present a theoretical analysis
of the excitation of a single trapped atom by HG modes. By expressing the HG
vector potentials in terms of their Bessel counterparts, written in the
paraxial approximation, we were able to derive simple analytical expressions
for the Rabi frequencies $\Omega^{\rm(HG)}$ of the photoinduced transitions.
These expressions allow one to analyze the dependence of $\Omega^{\rm(HG)}$
both on the parameters of HG beams and on the orientation of the external
magnetic field, used to define the quantization axis of a target ion. The
developed theory was applied to investigate the case of the
${}^{2}S_{1/2}\to\;^{2}F_{7/2}$ E3 transition in an 171Yb+ ion localized in
the beam’s center. Calculations, performed for experimentally feasible
parameters, helped us to determine optimal combinations of incident light
polarization and orientation of the magnetic field, resulting in rather large
Rabi frequencies. For HG10 and HG01 modes, we have estimated that
$\Omega^{\rm(HG)}$ lies in the range of tens of Hz, which can exceed
$\Omega^{(\text{LG})}$ of the previously used “standard” Laguerre-Gaussian
light for the same beam waist and power. This suggests that Hermite-Gaussian
beams can serve as a useful tool for studying highly-forbidden atomic
transitions.

In the present theoretical study we have used a “semiclassical” approach, Eq.
(15), to account for the delocalization of a target ion with respect to the
beam center. Even though this approach is expected to provide reliable
qualitative predictions for the dependence of the Rabi frequency on the
geometry of the experiment and on the state of the incident HG beam, it should
be probed against the fully quantum theory in which the center-of-mass motion
of a target atom is quantized. Similar studies have been performed for an
atomic target exposed to standing electromagnetic waves [13]. The fully
quantum analysis employing the time-dependent density matrix approach [34] is
currently underway and will be presented in a follow-up publication.

Following our theoretical analysis we propose a photonic integrated circuit
(PIC) frequency reference based on an octupole transition in Yb+ ions that can
be used in a PIC-based optical clock. Such a clock will be orders of magnitude
more compact compared to existing setups and have lower power consumption and
higher pointing stability, making it well suited for applications outside the
laboratory, including deep space navigation and geodetic measurements. As a
first experiment in this direction, we plan to use $x$-polarized HG10 or
$y$-polarized HG01 modes in a “clock-on-a-chip” setup. For lowest light
shifts, we plan to operate the clock with 173Yb+ which offers 10 times larger
transition matrix elements and thus a 100 fold reduction of the light shift on
its clock transitions [30, 35]. Since the linewidth of the atomic transition
is still narrower than the laser linewidth, depending on the laser noise
spectrum the suppression factor can be even larger.

## Acknowledgments

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – Project-ID 445408588 (SU 658/5-1) and under Germany’s
Excellence Strategy - EXC-2123 QuantumFrontiers - 390837967. This work was
also supported by the Max-Planck-RIKEN-PTB-Center for Time, Constants and
Fundamental Symmetries. We thank C.-F. Grimpe for the discussion about
integrated optics.

## References

  * [1] R. Heintzmann, M. G. L. Gustafsson, Nature Photon. 2009, 3, 362.
  * [2] M. Padgett, R. Bowman, Nature Photon. 2011, 5, 343.
  * [3] M. F. Andersen, C. Ryu, P. Cladé, V. Natarajan, A. Vaziri, K. Helmerson, W. D. Phillips, Phys. Rev. Lett. 2006, 97, 170406.
  * [4] J. Wang, J. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. Ren, Y. Yue, S. Dolinar, M. Tur, A. E. Willner, Nature Photon. 2012, 6, 488.
  * [5] M. Krenn, M. Huber, R. Fickler, R. Lapkiewicz, S. Ramelow, A. Zeilinger, Proc. Natl Acad. Sci. 2014, 111, 6243.
  * [6] H. Rubinsztein-Dunlop et al., J. Opt. 2017, 19, 013001.
  * [7] Y. Shen, X. Wang, Z. Xie, C. Min, X. Fu, Q. Liu, M. Gong, X. Yuan, Light Sci. Appl. 2019, 8, 90.
  * [8] C. T. Schmiegelow, J. Schulz, H. Kaufmann, T. Ruster, U. G. Poschinger, F. Schmidt-Kaler, Nat. Commun. 2016, 7, 12998.
  * [9] R. Lange, N. Huntemann, A. A. Peshkov, A. Surzhykov, E. Peik, Phys. Rev. Lett. 2022, 129, 253901.
  * [10] M. Solyanik-Gorgone, A. Afanasev, C. E. Carlson, C. T. Schmiegelow, F. Schmidt-Kaler, J. Opt. Soc. Am. B 2019, 36, 565.
  * [11] S. A.-L. Schulz, A. A. Peshkov, R. A. Müller, R. Lange, N. Huntemann, Chr. Tamm, E. Peik, A. Surzhykov, Phys. Rev. A 2020, 102, 012812.
  * [12] A. B. Mundt, A. Kreuter, C. Becher, D. Leibfried, J. Eschner, F. Schmidt-Kaler, R. Blatt, Phys. Rev. Lett. 2002, 89, 103001.
  * [13] A. R. Vasquez, C. Mordini, C. Vernière, M. Stadler, M. Malinowski, C. Zhang, D. Kienzler, K. K. Mehta, J. P. Home, Phys. Rev. Lett. 2023, 130, 133201.
  * [14] K. K. Mehta, C. D. Bruzewicz, R. McConnell, R. J. Ram, J. M. Sage, J. Chiaverini, Nature Nanotech. 2016, 11, 1066.
  * [15] K. K. Mehta, C. Zhang, M. Malinowski, T.-L. Nguyen, M. Stadler, J. P. Home, Nature 2020, 586, 533.
  * [16] R. J. Niffenegger, J. Stuart, C. Sorace-Agaskar, D. Kharas, S. Bramhavar, C. D. Bruzewicz, W. Loh, R. T. Maxson, R. McConnell, D. Reens, G. N. West, J. M. Sage, J. Chiaverini, Nature 2020, 586, 538.
  * [17] M. Babiker, D. L. Andrews, V. E. Lembessis, J. Opt. 2019, 21, 013001.
  * [18] B. A. Knyazev, V. G. Serbo, Phys.-Usp. 2018, 61, 449.
  * [19] A. Surzhykov, D. Seipt, V. G. Serbo, S. Fritzsche, Phys. Rev. A 2015, 91, 013403.
  * [20] M. E. Rose, _Elementary Theory of Angular Momentum_ , John Wiley & Sons, New York 1957.
  * [21] A. Afanasev, C. E. Carlson, C. T. Schmiegelow, J. Schulz, F. Schmidt-Kaler, M. Solyanik, New J. Phys. 2018, 20, 023032.
  * [22] B. Baghdasaryan, B. Böning, W. Paufler, S. Fritzsche, Phys. Rev. A 2019, 99, 023403.
  * [23] Edited by D. L. Andrews, M. Babiker, _The Angular Momentum of Light_ , Cambridge University Press, Cambridge 2013.
  * [24] G. J. Gbur, _Singular Optics_ , CRC Press, Boca Raton, FL 2017.
  * [25] M. Fox, _Quantum Optics: An Introduction_ , Oxford University Press, Oxford 2006.
  * [26] S. A. Akhmanov, S. Yu. Nikitin, _Physical Optics_ , Clarendon Press, Oxford 1997.
  * [27] V. A. Dzuba, V. V. Flambaum, M. S. Safronova, S. G. Porsev, T. Pruttivarasin, M. A. Hohensee, and H. Häffner, Nature Phys. 2016, 12, 465.
  * [28] C. Sanner, N. Huntemann, R. Lange, Chr. Tamm, E. Peik, M. S. Safronova, S. G. Porsev, Nature 2019, 567, 204.
  * [29] L. S. Dreissen, C.-H. Yeh, H. A. Fürst, K. C. Grensemann, T. E. Mehlstäubler, Nat. Commun. 2022, 13, 7314.
  * [30] R. Lange, A. A. Peshkov, N. Huntemann, Chr. Tamm, A. Surzhykov, E. Peik, Phys. Rev. Lett. 2021, 127, 213001.
  * [31] R. V. Popov, A. V. Maiorova, Opt. Spectrosc. 2017, 122, 366.
  * [32] M. Auzinsh, D. Budker, S. M. Rochester, _Optically Polarized Atoms: Understanding Light-Atom Interactions_ , Oxford University Press, Oxford 2010.
  * [33] W. M. Itano, J. Res. Natl. Inst. Stand. Technol. 2000, 105, 829.
  * [34] A. A. Peshkov, Y. M. Bidasyuk, R. Lange, N. Huntemann, E. Peik, A. Surzhykov, Phys. Rev. A 2023, 107, 023106.
  * [35] V. A. Dzuba, V. V. Flambaum, Phys. Rev. A 2016, 93, 052517.


# Data-driven detector signal characterization with constrained bottleneck
autoencoders

C.Jesús-Valls,11footnotetext: Corresponding author. T.Lux  F.Sánchez

###### Abstract

A common technique in high energy physics is to characterize the response of a
detector by means of models tuned to data which build parametric maps from the
physical parameters of the system to the expected signal of the detector. When
the underlying model is unknown it is difficult to apply this method, and
often, simplifying assumptions are made introducing modeling errors. In this
article, using a waveform toy model we present how deep learning in the form
of constrained bottleneck autoencoders can be used to learn the underlying
unknown detector response model directly from data. The results show that
excellent performance results can be achieved even when the signals are
significantly affected by random noise. The trained algorithm can be used
simultaneously to perform estimations on the physical parameters of the model,
simulate the detector response with high fidelity and to denoise detector
signals.

##  1 Introdution

In high energy physics (HEP), an usual analysis strategy is that of comparing
experimental measurements to synthetic datasets generated by means of a
simulator. Given a simulation model, the synthetic data can be tuned by
changing the underlying model parameters allowing to compute best fit points
and confidence intervals. However, to accurately model a system a precise
quantitative knowledge of its behavior is necessary. Sometimes, such a
knowledge is only partial, and the discrepancies between the simulated and the
measured datasets lead to systematic errors, decreasing the experimental
sensitivity.
In general, HEP simulations consist of two parts: The propagation of particles
in the detector and the detector response. Whereas the same modeling of
physics can be used in a large variety of experiments, the detector response
is often quite specific. Hence, dedicated characterization studies of a
detector are often necessary typically consisting of a multitude of smaller
characterizations tasks each describing the behavior of either a sub-detector
module or an individual channel. Modeling the detector signals is often non-
trivial due to the many effects that can arise in the detector response,
either in the signal production (e.g. quenching, de-excitations, etc), in the
signal transport (e.g. field non-uniformities, attenuation, diffusion,
reflections, spreading, etc) or in the sensor response (e.g. saturation,
after-pulsing, cross-talk, etc). Therefore, mismodeling often arises in the
detector simulation leading to systematic errors. To improve this, data-driven
solutions that characterize the signal with high fidelity are of great
interest as they can be used to replace, or improve, the inaccurate detector
response models. Due to this, a usual practice is that of adjusting laboratory
measurements with functions and use them to characterize the detector signal.
Although this is useful for a number of problems, this method is not quite
general, as it works only for scalar detector signals. If the signal is
instead an array, such as a waveform, the problem grows in complexity. Due to
this, simplifying assumptions must often be made which introduce modeling
errors. In this article, we present and discuss how constrained autoencoders
can be used to learn the underlying unknown detector response models directly
from data examples. We illustrate the potential of this method through the
study of a toy model producing signal waveforms and we discuss how the trained
algorithm can be used to perform parameter estimations, simulate the detector
response with high fidelity and to effectively denoise detector signals.

###  1.1 Problem definition

Consider the deterministic response of a detector as a transformation $f$ that
maps a set of physics parameters, $\vec{p}$, (e.g. track position, angle,
energy deposit, etc) into a set of detector signals $\vec{s}$, namely,

$f(\,\,\vec{p}\,\,)\mapsto\vec{s};\qquad
f^{-1}(\,\,\vec{s}\,\,)\mapsto\vec{p}\,.$ (1.1)

Then, our aim is to learn transformations $\hat{f}$ (and $\hat{f}^{-1}$) as
similar as possible to $f$ (and $f^{-1}$) directly from data with the only
requirement of knowing the set of $\vec{p}$ playing a role in the system and
being able to collect data samples of $\vec{s}$ for known $\vec{p}$
configurations. Noticeably, this single condition can often be satisfied, as
output signals $\vec{s}$ can be recorded in the laboratory while determining
$\vec{p}$ with the aid of additional instruments. An exemplifying toy model is
presented in the following subsection.

###  1.2 Toy problem

We will use a toy transformation $f$ given by:

$\displaystyle\vec{s}(x)$ $\displaystyle=f(\vec{p}=\\{\theta,\eta\\})(x)$
$\displaystyle=(1+0.5\theta)\mathcal{M}(-7|1+0.3\theta)+\mathcal{M}(7(\eta-1)|1.2+\theta)$
(1.2)

Where the shorthand notation $\mathcal{M}(\mu|c)$ has been used to define a
Moyal distribution [1] with location parameter $\mu$ and scale $c$, which
describes the energy loss of a charged relativistic particle due to ionization
of the medium and approximates the Landau distribution with a systematically
lower tail. This model takes two generic input physics parameters
$\vec{p}=\\{\theta,\eta\\}$. To give some physical intuition on the model it
has been inspired to resemble the response of a scintillating bar, similar to
the signals reported in Ref. [2]. Under this picture, the first (second) Moyal
distribution can be interpreted as accounting for the direct (reflected)
scintillating light reaching a photosensor in one end of the bar. Then, one
can think of the parameter $\theta$ as accounting for the track angle, related
to the light yield, and of $\eta$ as the distance at which the track has
crossed the bar with respect to the photosensor, introducing a time delay in
the measured signal. If $x$ is interpreted as time, a signal $s(x)$ is
measured in the detector at discrete consecutive instances of $x$. In this
way, the array of signals $\vec{s}$ constitutes a 1D waveform.
Now consider the following situation. The scintillating bar of the previous
example can be measured in the laboratory such that a collection of $\vec{s}$
can be recorded. For each of this signals, $\vec{p}$ is also determined with
the aid of additional instruments, e.g. additional position sensitive
detectors. However, the experimenter is unaware of the true functional form of
$f$, which can be in practice arbitrarily complex, but desires to learn $f$,
its inverse transformation, or both, out of the collected sets of $\vec{s}$
and $\vec{p}$.
In this article, we will use the toy model above to generate data instances of
$\vec{s}$ and $\vec{p}$ artificially, reproducing what could be collected in
the laboratory. Later, we will pretend that $f$ is unknown, and we will face
the problem of learning $\hat{f}$ and $\hat{f}^{-1}$ as close as possible to
$f$ and $f^{-1}$ directly from the synthetic data examples.

##  2 Autoencoders

Over the last years deep learning methods have become crucial for modern data
analysis and have started to play a important role in HEP. So far, most
solutions have been reached through supervised algorithms, particularly
applied to classification problems, with deep neural networks at the
forefront. Recently, deep generative models have aroused as a novel
unsupervised alternative, able to deal with both labeled and unlabeled
datasets and to learn sophisticated transformation between spaces of very
different dimensionality [3]. Autoencoders are a particular type of this novel
algorithms. In autoencoders (AEs) a real input of arbitrary dimensionality
$\mathbb{R}^{A}$ is mapped into a real space output of the same dimensionality
by following a set of transformations characterized by a bottleneck
architecture. The bottleneck, often referred to as the latent space, has a
reduced dimensionality $\mathbb{R}^{B}$, where typically $A\gg B$. As the
algorithm is trained to minimize the difference between the input and the
output distributions the autoencoder learns a compact representation of the
data which is encoded into the latent space. Due to this, the transformation
which maps the input into the latent space is named the encoder whereas the
transformation which maps the latent space into the output is named the
decoder.
Despite their recent appearance in the physics literature, AEs have already
been studied for several interesting HEP applications. On one hand, AEs can be
used to compress data by storing the much lighter latent space representation
[4, 5, 6]. On the other hand, AEs average out the noise from different data
examples during training and, therefore, they have great potential to denoise
signals [7, 8, 9]. AE transformations are greatly interesting by themselves.
If an AE is trained using non-anomalous signals, for instance produced with a
well understood simulation, the AE reconstruction error can later be used to
identify outliers and tag them as anomaly candidates. Thus, so far, anomaly
detection has been among the most widely explored applications involving AEs
in the HEP literature [10, 11, 12, 13, 14]. AEs can also be used to learn
signal or background distributions directly from data. In this applications,
AEs enable the generation of synthetic data by sampling in the latent space,
an option explored in several studies [15, 16, 17]. However, if the latent
space is unconstrained, no physical interpretation can be made of the latent
space. Due to this, architectures in where the bottleneck is constrained have
been studied in computer science, leading to the proposal of modified AEs,
such as the Bounded-Information-Bottleneck AE (BIB-AE) [18]. The main
difference with a standard AE is that in these modified algorithms additional
information in the form of labels is used in the training phase in order to
make the latent space interpretable. In HEP, this approach has been first
explored recently via an algorithm dubbed the end-to-end Sinkhorn AE. Its use
has been explored in the context of producing complex high fidelity simulation
outputs for calorimeters by replacing simulators with AEs trained in simulated
data [19, 20]. This has the main advantage of reducing the time and
computational cost of generating simulated events. In this article, we further
discuss the potential of latent space constrained autoencoders in HEP. In
particular, we focus on a powerful and unexplored application: its use for
data-driven detector signal characterization and, even more, how the same
trained algorithm can be used not only to generate signals but also to do
physical parameter estimations and signal denoising.

###  2.1 The constrained bottleneck autoencoder

We will use a modified AE, similar to the BIB-AE, with simple features which,
for generality, will be referred to as a constrained bottleneck autoencoder. A
sketch of the algorithm is presented in Figure 1. Its main features are as
follows: For training, events consisting of signals $\vec{s}$ associated to
physical configuration parameters $\vec{p}$ are given as the input to the
algorithm. The algorithm then performs a series of sequential transformations
which generate an output signal $\vec{s^{\prime}}$ with equal dimensionality
to $\vec{s}$. The bottleneck $\vec{p^{\prime}}$ corresponds to the centermost
transformation layer and has the same dimensionality as $\vec{p}$. With the
former settings, the network weights specifying the encoder and decoder
transformations $\hat{f}^{-1}$ and $\hat{f}$, corresponding to our target
interest (see Sec 1.1), can be optimized by minimizing the loss over the set
of input examples. Such loss consists of the sum of two terms. On one hand,
the mean squared error (MSE) between $\vec{s}$ and $\vec{s^{\prime}}$ which
accounts for the similarity of the input and output signals. On the other
hand, the MSE between $\vec{p}$ and $\vec{p^{\prime}}$ which accounts for the
similarity between the latent space and the physical configuration parameters.

Figure 1:  Sketch of the constrained bottleneck autoencoder under discussion.

###  2.2 Architecture and implementation

The former architecture design has been driven by simple principles. In the
first place, using convolutional transformation allows to identify patterns
arising in the correlations among neighbor signal values. The number of
iterations, filters, kernel size and stride were chosen after a set of
dedicated architectural tests on the overall performance of the network. In
the second place, the use of fully connected layers allows to combine the
information from the different convolution filters and provides to the
algorithm enough flexibility in the encoder and decoder transformations. Such
flexibility is crucial, as our AE has a constrained bottleneck, and therefore,
learning $\hat{f}$ and $\hat{f}^{-1}$ is significantly more challenging than
in the unconstrained case. In this sense, stacking N-times the former layers
allows to organically increase the complexity of the network, either by
increasing the number of filters, the number of repetitions, or both.

##  3 Methodology

###  3.1 Dataset generation

A dataset has been generated by producing a set of 50k statistically
independent $\vec{s}$. Each instance of $\vec{s}$ has been generated using
uncorrelated flat random values of $\theta\in[0,1]$ and $\eta\in[0,1]$ and by
evaluating $f$, as defined in Eq. 1.2, a total of $2^{10}=1024$ equidistant
times along an interval of $x\in[-10,10]$. In this way, the signals $\vec{s}$
are waveforms of 1024 values. For training purposes the waveforms have been
normalized. This can be achieved in general by applying the change of
variable: $\vec{s}=(\vec{s}+\textup{offset})/\textup{scale}$. For convenience,
we have used offset=0.1 and scale=0.8 such that after adding noise to the
waveforms (explained below), the signals are made of values largely contained
in the interval [0,1]. An illustration of the normalized waveforms provided by
the toy model in Eq. 1.2 is presented in Figure 2.

Figure 2:  Functional normalized shape of $f(\vec{p})$ as described in Eq 1.2.
A colormap of $\theta$ for a series of slices of $\eta$ is used to illustrate
the dependency on both parameters.

In a realistic case, the waveforms would be affected by random fluctuations.
To show the robustness of the method to random noise and to illustrate its
denoising potential prior to the signal normalization each value of $\vec{s}$
has been varied independently, adding to it the result of random sampling from
a Gaussian with the mean at the original signal value and a sigma of 0.05. At
the normalization scale=0.8, this level of noise corresponds to $6.25\%$ of
the maximum possible signal amplitude. Example results of a set of 15 randomly
chosen instances of $\vec{s}$ compared to their associated $f(\vec{p})$ are
presented in Figure 3.

Figure 3:  Examples of 15 random noisy waveforms in the test dataset compared
to its original functional shape, specified by Eq. 1.2.

###  3.2 Training

For training and testing the system the dataset of 50k examples has been split
into two separate datasets: A training dataset with 40k signals and a test
dataset with the remaining 10k events. The autoencoder has been trained for 10
epochs (complete cycles to all the data), using the Adam optimizer [21] with
an initial learning rate of $5\cdot 10^{-4}$ and a mini-batch size of 32
(which determines the amount of examples that are processed before updating
the network parameters). During training the network weights are optimized by
minimizing the loss function, which as earlier anticipated in Figure 1,
consists on the addition of the minimum squared error of the signal
reconstruction accuracy and the latent space matching accuracy. The training
for which the results are reported took about 50 min in a standard
computer2222.3 GHz dual-core Intel Core i5..

##  4 Results and discussion

The performance of the trained algorithm can be evaluated for different
purposes. In the first hand, the quality of the learned
$\hat{f}^{-1}(\vec{s})\mapsto\vec{p^{\prime}}$ transformation is relevant to
do estimates of the physical parameters $\vec{p}$ using signals as input. On
the second hand the quality of the learned
$\hat{f}(\vec{p})\mapsto\vec{s^{\prime}}$ is relevant to model the detector
response and to produce realistic signals $\vec{s}$, e.g. if the algorithm is
embedded in a Monte Carlo simulation. Notably, since the algorithm learns from
multiple examples the random noise is averaged out, such that
$\hat{f}(\vec{p})$ generates noiseless signals $\vec{s^{\prime}}$, even if it
learns from noisy examples $\vec{s}$. Due to this, the transformation
$\hat{f}^{-1}(\hat{f}(\vec{s}))\mapsto\vec{s^{\prime}}$ can be used to denoise
signals. We study each of these tasks separately and discuss them below.

###  4.1 Reconstruction of the physical parameters

Figure 4:  Distribution of the two reconstructed physical parameters
($\theta^{\prime}$ and $\eta^{\prime}$) of the toy model for the whole test
dataset. The reconstructed values are shown in slices of true $\eta$ and
colored with the true value of $\theta$, allowing to visualize the latent
space. The horizontal black lines show the constrains on the true values of
$\eta$.

Figure 5:  2D distribution of the reconstruction error for $\theta$ and $\eta$
and its 1D marginal error distributions including Gaussian fits. The
correlation of the two errors is also presented.

The distribution of the reconstructed physical parameters,
$\vec{p}=\\{\theta,\eta\\}$, associated to the signals for the whole test
dataset of 10k events is Figure 5. The results are presented in slices of true
$\eta$ and in a colormap of $\theta$ allowing to qualitatively understand the
2D latent space distributions in the whole phase-space. As expected, the
encoder learns a transformation that maps signals $\vec{s}$ associated to
transformations determined by values of $\vec{p}=\\{\theta,\eta\\}$ to a
latent space $\vec{p^{\prime}}=\\{\theta^{\prime},\eta^{\prime}\\}$ where the
typical distance between $\vec{p}$ and $\vec{p}^{\prime}$ is small. This is
explicitly shown in Figure 5 where the errors in each of the reconstructed
physical parameters is presented as a 2D plot showing pairs of reconstructed
$\\{\theta^{\prime},\eta^{\prime}\\}$ and their 1D marginal distributions.
Notably, the errors are Gaussian distributed, and show negligible bias and a
very narrow standard deviation. The width of the distribution is different for
the two parameters. This is expected, as the two parameters affect differently
the waveform shape. In particular, looking into Figure 2 it is clear that a
wide range of $\theta$ can render very similar waveforms, as for instance for
small $\eta$ values. From these results on the reconstructed latent space two
major conclusions can be extracted. In the first place the algorithm is able
to effectively learn the transformation
$\hat{f}^{-1}(\vec{s})\mapsto\vec{p^{\prime}}$, allowing to do parameter
estimations directly from detector signals. In the second place, no outliers
are observed and the error distribution is Guassian allowing to easily and
reliable associate errors to the parameters estimated by means of this method.

###  4.2 Signal generation

Figure 6:  Functional shape of $\hat{f}(\vec{p})$ as learned by the algorithm.
A colormap of $\theta$ for a series of slices of $\eta$ is used to illustrate
the dependency on both parameters. Figure 7:  Subtraction of Figure 2 and
Figure 6 to illustrate the difference of $f(\vec{p})-\hat{f}(\vec{p})$.

The trained model can be used to characterize the transformation across the
whole phase space, as presented in Figure 6. In the ideal case, this scan
should be equivalent to that earlier presented in Figure 2. Eye inspection
allows to draw qualitative conclusions. In the first place, overall, the
algorithm has learned a transformation $\hat{f}(\vec{p})$ directly from data
examples which generally resembles the true one, and the predicted signals do
not show worrisome artifacts or instabilities. This is remarkable given that
the transformation is learned from noisy data examples but the algorithm
directly provides a denoised output. To draw quantitative conclusions, in
Figure 7 we present a plot of the difference of the true and the learned
model, namely a plot of $f(\vec{p})-\hat{f}(\vec{p})$. For the majority of the
phase-space the variations are very small, typically well contained within 2%
of the maximum amplitude. The largest differences are below 4% of the maximum
amplitude and are therefore sensibly smaller than the simulated noise level.

###  4.3 Signal denoising

Figure 8:  Left panel: Distributions of the signal reconstruction error
including various quantiles for one every sixteen waveform positions for the
whole test dataset of 10000 waveforms. Right panel: Aggregated distribution of
all 1024$\times$10000 signal reconstruction errors (cyan), and its 1024
medians (orange). The signal reconstruction error is defined as the difference
between the normalized noiseless signals provided by Eq. 1.2 and the signals
generated by the transformation
$\hat{f}^{-1}(\hat{f}(\vec{s}))\mapsto\vec{s^{\prime}}$. Figure 9:  The input
noisy signal $\vec{s}$ is compared to the denoised generated signal
$\vec{s^{\prime}}$ and to the true shape of the toy model $f(\vec{p})$ for a
total of 35 examples in the test dataset. The examples are chosen by their
mean squared reconstruction error percentile. Percentile 0 (100) corresponds
to the best (worse) reconstructed signal in the whole test dataset of 10k
waveforms.

Finally, as we have learned both $\hat{f}(\vec{s})\mapsto\vec{p^{\prime}}$ and
$\hat{f}^{-1}(\vec{p})\mapsto\vec{s^{\prime}}$ we can study the reconstruction
of $\hat{f}^{-1}(\hat{f}(\vec{s}))\mapsto\vec{s^{\prime}}$. This has two main
advantages. On one hand the results of the composite transformation map noisy
signals to denoised signals. In the second hand, the errors in the
reconstruction of physical parameters and in the generation of signals are
correlated, and in consequence the composite transformation might have a
superior performance than the two individual transformations.
To evaluate the quality of the reconstructed denoised signals
$\vec{s^{\prime}}$ from noisy signals $\vec{s}$, the output of
$\hat{f}^{-1}(\hat{f}(\vec{s}))\mapsto\vec{s^{\prime}}$ has been compared to
$f(\vec{p})$ for the 10k events in the test dataset. The results are presented
in Figure 8. In the left of the Figure, the signal reconstruction error
distributions for all waveform positions multiple of 16 are presented. The
decision of showing only one in every 16 positions is driven by a clearer
visualization, reducing the number of distributions from 1024 to 64. To show
information for all the waveform positions, in the right part of the Figure
aggregated results for all the 1024 positions are shown. These aggregated
results consists of a distribution of the 1024 median positions in each
waveform position (orange), and the 1024$\times$104 values reconstructed in
any position for all events (cyan).
The results allow to draw numerous conclusions. In the first place, across all
waveform positions the errors in the median are remarkably small, with a
typical bias (standard deviation) similar to 0.3% (0.2%) of the maximum
amplitude. The error in all reconstructed signals, is slightly wider but still
of excellent quality with a typical bias of (standard deviation) 0.3% (0.7%).
In the second place the reconstructed signals in all positions are densely
packed around around the median without any significant outlier. The maximum
error for the 64 positions under consideration is comparable with noise
fluctuations and a very smooth and well controlled signal prediction by the
network across the whole test dataset is observed.

To illustrate the accuracy and signal fidelity of the results presented in
Figure 8 all the 10k events in the test dataset have been sorted by their mean
squared signal reconstruction error. Waveforms at different error percentiles
are presented in Figure 9, being the percentile 0 (percentile 100) the event
with minimum (maximum) signal reconstruction error. The results show excellent
prediction fidelity across all percentiles with very minor differences even
for percentile 100.

##  5 Conclusions

An algorithm consisting of a constrained bottleneck autoencoder has been
presented as a reliable solution for the task of characterizing
multidimensional signals directly from data. In particular, a toy problem has
been presented and used to generate data. This data has been distorted by
random noise and used to train and test the algorithm, blind to the toy model
underlying functional form. Overall, the results show excellent performance in
learning the true distribution mapping physical parameters into signals and
the inverse of this distribution.
We have shown how to use the trained algorithm to do parameter estimations
directly from signals and how to associate errors to this estimates. We have
presented how, given some input physical parameters, the algorithm can be used
to model signals with good precision and hence it is a suitable candidate to
characterize the detector response in Monte Carlo simulations. Finally, we
have shown the excellent potential of this algorithm in denoising signals.

Each of the former results open potential applications in HEP were the
treatment and analysis of multidimensional signals, such as waveforms, is
widely spread. In addition, the excellent accuracy this method in mapping
signals to signals allows to envision compression applications in which the
much lighter bottleneck representation is stored instead of the full signal.

##  6 Acknowledgments

C. Jesús-Valls acknowledges fruitful discussions with S.Alonso-Monsalve,
C.Grimal-Bosch, S. Bordoni and L. Uboldi. C. Jesús-Valls and T. Lux
acknowledge funding from the Spanish Ministerio de Economía y Competitividad
(SEIDI-MINECO) under Grants No. PID2019-107564GB-I00 and SEV-2016-0588. IFAE
is partially funded by the CERCA program of the Generalitat de Catalunya. F.
Sánchez acknowledge the Swiss National Foundation Grant No. 200021_85012.

## References

  * [1] J. Moyal, “Xxx. theory of ionization fluctuations,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, vol. 46, no. 374, pp. 263–280, 1955.
  * [2] A. Korzenev et al., “A $4\pi$ time-of-flight detector for the ND280/T2K upgrade,” JINST, vol. 17, no. 01, p. P01016, 2022.
  * [3] L. Ruthotto and E. Haber, “An introduction to deep generative modeling,” GAMM-Mitteilungen, p. e202100008, 2021.
  * [4] G. Di Guglielmo et al., “A reconfigurable neural network ASIC for detector front-end data compression at the HL-LHC,” IEEE Trans. Nucl. Sci., vol. 68, p. 2179, 2021.
  * [5] Y. Huang, Y. Ren, S. Yoo, and J. Huang, “Efficient data compression for 3d sparse tpc via bicephalous convolutional autoencoder,” in 2021 20th IEEE International Conference on Machine Learning and Applications (ICMLA), pp. 1094–1099, IEEE, 2021.
  * [6] B. Yoon, C. C. Chang, G. T. Kenyon, N. T. T. Nguyen, and E. Rrapaj, “Prediction and compression of lattice QCD data using machine learning algorithms on quantum annealer,” PoS, vol. LATTICE2021, p. 143, 2021.
  * [7] P. Ai, D. Wang, G. Huang, N. Fang, D. Xu, and F. Zhang, “Timing and characterization of shaped pulses with MHz ADCs in a detector system: a comparative study and deep learning approach,” JINST, vol. 14, no. 03, p. P03002, 2019.  [Erratum: JINST 14, E03001 (2019)].
  * [8] H. Shen, D. George, E. A. Huerta, and Z. Zhao, “Denoising gravitational waves with enhanced deep recurrent denoising auto-encoders,” in ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3237–3241, 2019.
  * [9] M. Erdmann, F. Schlüter, and R. Smida, “Classification and Recovery of Radio Signals from Cosmic Ray Induced Air Showers with Deep Learning,” JINST, vol. 14, no. 04, p. P04005, 2019.
  * [10] M. Farina, Y. Nakai, and D. Shih, “Searching for New Physics with Deep Autoencoders,” Phys. Rev. D, vol. 101, no. 7, p. 075021, 2020.
  * [11] J. Hajer, Y.-Y. Li, T. Liu, and H. Wang, “Novelty Detection Meets Collider Physics,” Phys. Rev. D, vol. 101, no. 7, p. 076015, 2020.
  * [12] A. Blance, M. Spannowsky, and P. Waite, “Adversarially-trained autoencoders for robust unsupervised new physics searches,” JHEP, vol. 10, p. 047, 2019\.
  * [13] M. Crispim Romão, N. F. Castro, and R. Pedro, “Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders,” Eur. Phys. J. C, vol. 81, no. 1, p. 27, 2021.  [Erratum: Eur.Phys.J.C 81, 1020 (2021)].
  * [14] T. Cheng, J.-F. Arguin, J. Leissner-Martin, J. Pilette, and T. Golling, “Variational Autoencoders for Anomalous Jet Tagging,” arXiv:2007.01850, 7 2020.
  * [15] P. Holl, L. Hauertmann, B. Majorovits, O. Schulz, M. Schuster, and A. J. Zsigmond, “Deep learning based pulse shape discrimination for germanium detectors,” Eur. Phys. J. C, vol. 79, no. 6, p. 450, 2019.
  * [16] C.-H. Liao and F.-L. Lin, “Deep generative models of gravitational waveforms via conditional autoencoder,” Phys. Rev. D, vol. 103, no. 12, p. 124051, 2021.
  * [17] J. W. Monk, “Deep Learning as a Parton Shower,” JHEP, vol. 12, p. 021, 2018.
  * [18] S. Voloshynovskiy, M. Kondah, S. Rezaeifar, O. Taran, T. Holotyak, and D. J. Rezende, “Information bottleneck through variational glasses,” arXiv preprint arXiv:1912.00830, 2019.
  * [19] E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol, and K. Krüger, “Getting High: High Fidelity Simulation of High Granularity Calorimeters with High Speed,” Comput. Softw. Big Sci., vol. 5, no. 1, p. 13, 2021.
  * [20] E. Buhmann et al., “Fast and Accurate Electromagnetic and Hadronic Showers from Generative Models,” EPJ Web Conf., vol. 251, p. 03049, 2021\.
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (Y. Bengio and Y. LeCun, eds.), 2015.