# Doppler Compensated Cavity For Atom Interferometry

Rustin Nourshargh1, Sam Hedges1, Mehdi Langlois1, Kai Bongs1, Michael
Holynski1∗

1MUARC, School of Physics and Astronomy, University of Birmingham
Edgbaston, Birmingham, B15 2TT, UK
∗Corresponding author. E-mail<EMAIL_ADDRESS>
> We propose and demonstrate a scheme to enable Doppler compensation within
> optical cavities for atom interferometry at significantly increased mode
> diameters. This has the potential to overcome the primary limitations in
> cavity enhancement for atom interferometry, circumventing the cavity
> linewidth limit and enabling mode filtering, power enhancement, and a large
> beam diameter simultaneously. This approach combines a magnified linear
> cavity with an intracavity Pockels cell. The Pockels cell introduces a
> voltage tunable birefringence allowing the cavity mode frequencies to track
> the Raman lasers as they scan to compensate for gravitationally induced
> Doppler shifts, removing the dominant limitation of current cavity enhanced
> systems. A cavity is built to this geometry and shown to simultaneously
> realize the capability required for Doppler compensation, with a 5.04 mm
> $1/e^{2}$ diameter beam waist and an enhancement factor of $>$5x at a
> finesse of 35. Furthermore, this has a tunable Gouy phase, allowing the
> suppression of higher order spatial modes and the avoidance of regions of
> instability. This approach can therefore enable enhanced contrast and longer
> atom interferometry times while also enabling the key features of cavity
> enhanced atom interferometry, power enhancement and the reduction of
> aberrations. This is relevant to future reductions in the optical power
> requirement of quantum technology, or in providing enhanced performance for
> atom interferometers targeting fundamental science.

## Introduction

Light-pulse atom interferometry is a technique employing atomic superposition
states as a sensitive probe. Laser pulses act as atomic beam splitters,
transferring atoms into a superposition of states where they are left to
evolve. Further pulses recombine the superposition, and the final atomic state
populations serve as an exquisite record of the environment in which they
evolved. The samples, clouds of laser cooled atoms, are interrogated to
measure quantities such as gravity, accelerations, and rotations (?, ?, ?, ?,
?). The sensitivity of these instruments is maximized when the interrogation
pulses are completed with high fidelity. This requires laser beams with a
large diameter compared to the cloud size, high intensities, and aberration-
free optical wavefronts with the latter being a leading cause of uncertainty
(?, ?, ?, ?, ?). Improvements in instrument sensitivity are achieved by
increasing the interferometer space-time area. This can be realized through
longer interferometry times or increasing the momentum transferred to the
atoms. Interferometry time is typically limited by the available free-fall
distance (?), and the thermal expansion of the cloud. Momentum transfer is
limited by the atom-optic beam splitter fidelity which depends on the
homogeneity of the beam profile, and the pulse parameters which are in
practice limited by the available laser power (?). Using an optical cavity to
provide resonant power enhancement and spatial mode filtering of the atom
interferometry beams offers considerable promise in addressing these
challenges. Existing cavity enhanced atom interferometers have demonstrated
the exciting potential of the technique (?, ?). Challenges implementing large
cavity modes in a geometrically stable way and Doppler shift limited
interrogation times have so far made it difficult to deliver on all of the
advantages. Overcoming these challenges would allow cavity enhancement to
fully benefit atom interferometers, and have the potential to provide
significant advantages to a wide range of experiments, including within
fundamental physics (?, ?, ?, ?, ?, ?, ?, ?, ?) and practical applications (?,
?, ?, ?).

During the atom interferometry sequence changes in the atomic cloud
velocities, due to gravity or photon recoils, result in a Doppler shift of the
laser frequencies as seen by the atoms. Accurate control of the interferometer
phase requires compensation of this shift by chirping the laser frequency. In
this work we consider an atom interferometer employing two-photon Raman
transitions on the D2 line of 87Rb, requiring the frequency difference between
the two Raman lasers to be chirped at
$\alpha_{chirp}=25.1~{}\mathrm{MHz}\,\mathrm{s}^{-1}$. However, operation
within an optical cavity places restrictions on the laser frequencies, with
implications for atom interferometry as outlined in Materials and Methods.
Whilst multiple input frequencies separated by the cavity free spectral range
$\nu_{FSR}$ (or multiples thereof) can be resonant simultaneously, standard
cavity manipulations do not allow the spacing between resonant modes to be
adjusted dynamically. As a result the frequencies cannot be simultaneously
resonant with the cavity and compensate the Doppler shift at significant free-
fall times, limiting the achievable space-time area for cavity enhanced atom
interferometry. This has placed limits on the performance of existing cavity
enhanced interferometers (?, ?, ?). An approach for addressing this has been
recently demonstrated (?).

In addition, cavity enhanced atom interferometry with conventional two mirror
cavities has been limited by small mode diameters, comparable to the
millimeter scale of the atom cloud. This results in dephasing as atoms in the
cloud sample different Rabi frequencies at different locations across the
cavity mode. This causes a corresponding loss of contrast in the atom
interferometer output. In free space atom interferometers a beam to cloud size
ratio greater than ten is generally used. Realizing large modes in an optical
cavity is a challenging but tractable problem. Recent advancements have
yielded a marginally stable cavity with a mode diameter of 1 cm (?) but
marginal stability causes higher order mode degeneracy preventing effective
spatial mode filtering, one of the key benefits sought through cavity
enhancement.

In this work we demonstrate a scheme for enabling Doppler compensation while
achieving large mode diameters and the ability to operate in a stable regime
to provide spatial mode filtering. This is achieved through the combination of
a magnified linear cavity and an intracavity Pockels cell. This is the first
example of a scheme for simultaneously overcoming two of the main limitations
in cavity enhanced atom interferometry, limited free-fall times and small mode
volumes.

Figure 1: (A): Schematic overview of the cavity, not to scale. The large mode
interaction region is the length $L_{3}$. The small waist and in-coupling mode
is on mirror $m_{1}$ which is planar, as is $m_{2}$. The Pockels cell is
within $L_{1}$. The two lenses have focal lengths $f_{1}=12~{}\mathrm{mm}$ and
$f_{2}=150~{}\mathrm{mm}$ yielding 12.5$\times$ magnification. (B): atoms
(purple) experience Doppler shifted laser frequencies when moving with respect
to the cavity. (C): The linearly polarized Raman beams, $R_{1}$ and $R_{2}$
are aligned with the orthogonal principle axes of the Pockels cell, $n_{x}$
and $n_{y}$ respectively. (D): The Pockels effect enabling tracking of a 40
MHz frequency shift in a simulated cavity with a finesse of 35 and $\nu_{FSR}$
of 217 MHz. $P_{1}$ and $P_{2}$ are the resonant frequencies of the two
polarizations within the cavity, which track the chirping Raman lasers $R_{1}$
and $R_{2}$ in frequency as they compensate for the freefall Doppler shift. U
is the uncompensated resonant frequency.

## Results

Here the magnified geometry and intracavity Pockels cell approach is described
and demonstrated. A $12.5\times$ magnification intracavity telescope produces
a Gaussian beam with a $1/e^{2}$ diameter of $5.04~{}\mathrm{mm}$ in a
geometrically stable configuration whilst also enabling Gouy phase tuning to
maximize higher order mode suppression. An intracavity Pockels cell induces
voltage controlled birefringence, lifting the degeneracy of orthogonally
polarized cavity modes and enabling independent tuning of their resonant
frequencies. Agile control of the Pockels cell birefringence allows the cavity
resonances to track the interferometry lasers in frequency as they chirp to
compensate for Doppler shifts.

### Mode diameter

Producing large mode diameters within a cavity is technically challenging, and
generally requires either a long cavity, or results in higher order mode
degeneracy and marginal stability. We avoid these problems with an intracavity
telescope focused near infinity in a two mirror cavity, as shown in figure 1.
Planar mirrors were selected due to wider availability, leaving five free
parameters which must be optimized. Of these $L_{2}$ and $L_{1}$ are
particularly sensitive. An ABCD-matrix formalism (?) is used to evaluate the
mode radius on mirror-2 ($m_{2}$) over this parameter space with results shown
in figure 2. Optimized parameters are found to be $L_{3}=390~{}\mathrm{mm}$
and focal lengths $f_{1}=12~{}\mathrm{mm}$ and $f_{2}=150~{}\mathrm{mm}$.

We achieve a beam waist of $2w=5.04~{}\mathrm{mm}$ on $m_{2}$, as measured
with a CCD. Around the central set of lens spacings ($L_{2}$), this waist is
remarkably insensitive to perturbations in the telescope length. A beam
diameter fluctuation of 1% would require a change in $L_{2}$ of more than
$100~{}\mu\mathrm{m}$. The intrinsic stability of this cavity geometry is such
that we were able to construct it from standard optomechanics mounted to an
aluminum rail.

Figure 2: The mode radius in the large diameter section of the cavity is shown
for different values of $L_{1}$ and $L_{2}$. At the edges of the surface the
radius diverges as the cavity reaches the limits of stability ($m\to\pm 1$).
Specific slices through the surface illustrate the insensitivity of the mode
radius to perturbations in $L_{2}$. In conventional designs with large modes
the behavior is similar to that at the edges of this surface. Owing to the
large beam diameter and correspondingly long Rayleigh range, 5.04 mm and 25.6
m respectively, these results are almost completely independent of $L_{3}$
which is fixed at $L_{3}=0.39$ m. Contours representing the higher order mode
degenerate values of Gouy phase are shown in white, and should be avoided to
improve spatial mode performance. The cavity described in this work is
indicated by the red spot on each plot. The plot was generated by calculating
the ABCD matrices for $10^{6}$ different cavities to select optimal cavity
parameters.

### Spatial mode filtering

An optical cavity only supports spatial modes with well defined wavefront
profiles. By coupling into the fundamental mode of the cavity we are able to
exploit this property to filter unwanted spatial modes and achieve flat
optical wavefronts. To provide such filtering effectively, the cavity design
must avoid frequency degeneracy between higher order modes and the fundamental
mode. If this is not avoided parasitic coupling to these higher orders will
produce significant wavefront distortion. This can be resolved through tuning
the Gouy phase. Higher order spatial modes acquire Gouy phase more rapidly
than the fundamental Gaussian mode and have different resonant frequencies.
The magnitude of the Gouy phase acquired over a round trip is given by

$\phi_{p,q}=(p+q+1)\cos^{-1}(m)$ (1)

for a higher order mode parameterized by indices $p,q$ (e.g. $TEM_{p,q}$) in
terms of the half-trace of the ABCD matrix describing the cavity, $m$ (?).

Higher order mode degeneracy will occur when $\phi_{p,q}$ is close to a
multiple of $\pi$. Adjusting the spacing of the telescope allows the half-
trace $m$, and hence the Gouy phase, to be tuned freely over the full range.
The Gouy phase is tuned to a value reducing higher order mode degeneracy,
whilst maintaining only moderate sensitivity to perturbations in length.
Cavity parameters which result in degeneracy with the first five higher order
modes are shown in figure 2 (three contours and the two boundaries of the
surface). We select $m=0.31$ resulting in $\phi=1.25$ avoiding the degenerate
values of $\pi/2$ and $\pi/3$. Equivalent Gouy phase adjustments are possible
if the lenses are replaced by curved mirrors, reducing round trip loss but
introducing astigmatism (?).

### Doppler compensation using the Pockels Effect

Maintaining cavity resonance whilst independently chirping two input laser
frequencies is achieved through the use of an intracavity electro-optic
crystal. An intracavity Pockels cell induces a voltage controlled
birefringence in the cavity, increasing the optical path length in the slow
axis with respect to the fast axis. This lifts the degeneracy between the two
orthogonally polarized cavity modes and provides a differential shift of their
resonant frequencies allowing them to be independently tuned, see figure 3.
Orthogonal linearly polarized Raman beams implement the atom optics, and the
polarization axes are aligned with the principle axes of the Pockels cell, see
figure 1C. The frequency shift is linear in applied voltage. This induces a
$\lambda/2$ shift in the resonance conditions of $\nu_{FSR}/2$ in a single
pass traveling wave cavity and $\nu_{FSR}$ in a double pass standing-wave
cavity as presented in this work. The slew rate for Doppler compensation is
given by

$\frac{dV}{dt}=\frac{2}{N_{passes}}\frac{\alpha_{chirp}}{\nu_{FSR}}V_{\lambda/2}$
(2)

where the half-wave voltage depends on the Pockels cell crystal,
$V_{\lambda/2}=3.3~{}\mathrm{kV}$ for our system, and in the standing wave
geometry employed here $N_{passes}$ is two. The required slew rate for Doppler
compensation of rubidium in free-fall
($\alpha_{chirp}=25.1~{}\mathrm{MHz~{}s}^{-1}$), is therefore easily
achievable.

Pockels cells are bi-polar, meaning that sweeping the voltage from
$-V_{\lambda/2}\to V_{\lambda/2}$ allows the relative frequency to be varied
by $\nu_{\mathrm{FSR}}$ in a traveling wave cavity and $2\nu_{\mathrm{FSR}}$
in a standing wave cavity. If longer sweeps are required, the resonance
condition can be maintained by reversing the polarity of the applied voltage
between pulses and noting that the longitudinal modes now in use are separated
by an additional free spectral range (in a traveling wave configuration a
double length Pockels cell is required to achieve arbitrary frequency shifts).

Figure 3: Increasing the voltage applied to the Pockels cell induces
birefringence into the cavity. The orthogonally polarized cavity modes are
separated in frequency, from the degenerate value at -160 V to 114 MHz at 1600
V applied. The initial offset compensates for residual birefringence in the
crystal. Ramping the voltage at $390~{}\mathrm{V\,s^{-1}}$ compensates for the
$25.1~{}\mathrm{MHz\,s^{-1}}$ Doppler shift.

An intracavity Pockels cell also relaxes the restrictions on cavity lengths.
In previous work, fine tuning of the free spectral range $\nu_{FSR}$ was
required to ensure both Raman beams (separated by $\nu_{HFS}$) are resonant
with the cavity simultaneously (?). This is satisfied when
$\nu_{HFS}=N\nu_{FSR}$, that is, the frequency splitting between Raman beams
is an integer multiple of $\nu_{FSR}$. By allowing the initial Pockels cell
voltage to take a non-zero value, resonance can be achieved for both Raman
beams without restricting $\nu_{FSR}$. In a system where one of the Raman
frequencies is generated by electro-optic modulation, having a non-integer
$\nu_{FSR}$ separation carries the additional benefit of suppressing any
additional frequencies present during the atom interferometry sequence. For
instance, for electro-optically generated Raman beams this would suppress
unwanted sidebands which would otherwise cause AC-stark shifts and parasitic
transitions. Previously such sidebands would have been co-resonant. The degree
of suppression will vary as the Pockels cell is scanned and careful choice of
pulse timing, cavity $\nu_{FSR}$, and Gouy phase can be used to allow these
co-resonances to be avoided.

This technique compensates for Doppler shifts and allows the use of
arbitrarily large cavities and at high finesse. The photon lifetime in this
cavity is $\sim 50~{}\mathrm{ns}$, which for typical pulse durations of
$1-10~{}\mathrm{\mu s}$ will not cause significant temporal pulse distortion.
For longer, high finesse cavities with a photon lifetime comparable with the
pulse duration, distortion of the temporal profile of the pulses becomes
significant and further constraints apply (?, ?).

### Power enhancement

Resonant power enhancement is one of the key benefits we realized with the
cavity. This proof of principle system demonstrated a power enhancement of
over 5x with a finesse of $\mathcal{F}=35$, using standard components. This
has been achieved despite the cavity mode encountering 26 optical surfaces and
passing through intracavity optics. Whilst modest for a typical cavity, this
is a dramatic enhancement in the context of atom interferometry which could
further be improved through optimisation of the optical scheme and components.

## Discussion

We demonstrate a scheme for cavity enhanced atom interferometry that for the
first time combines a large mode diameter with the capability to provide
Doppler compensation, reduction of wavefront aberration through spatial mode
filtering and power enhancement. This overcomes the two key barriers
preventing high performance use of cavity enhanced atom interferometry, while
realising the two primary benefits. This is achieved while also being robust
to changes in cavity length, avoiding the extreme dimensional tolerances
encountered at the edges of geometric stability. The techniques described can
enable large-mode, high-finesse cavity experiments capable of compensating for
arbitrarily large Doppler shifts due to long free-fall times, and the
parameters required to enable large momentum transfer orders. Furthermore,
this scheme could be applied to state-of-the-art long baseline atom
interferometers, with 10 m baselines and 2 s freefall times (?), while still
remaining within the limitations imposed by cavity lifetime elongation (?). As
such, this scheme is highly applicable to the full range of applications for
atom interferometry. For fundamental physics, the scheme could be applied to
realize increased optical intensities and reductions in wave-front aberration
to enhance space-time area through large momentum transfer without imposing
limitations on free-fall time. For quantum technology sensors the scheme could
provide reductions in the required laser power, leading to cheaper sensors or
operation on power constrained platforms such as satellites, while also
providing a route to increased sensitivity through robust cavity enhancement
for large momentum transfer.

## Materials and Methods

### Cavity linewidth limit

Atoms in an atom interferometer aligned with gravity experience Doppler
shifted laser frequencies during freefall. Chirping the lasers to compensate
for the resulting Doppler shift has the undesired effect of moving both laser
frequencies away from cavity resonance. This cavity linewidth limit results in
reduced circulating intensity and prevents beam splitters entirely at long T
times. We arbitrarily define a maximum interrogation time $T_{max}$, to be
that at which the Doppler shift has reduced the circulating intensity to half
of it’s resonant value for the final beam splitter pulse ($2T_{max}$ total
free fall).

$T_{max}=\frac{1}{2}\frac{\nu_{FSR}}{\alpha_{chirp}\mathcal{F}}=\frac{1}{2}\frac{c}{\alpha_{chirp}L_{RT}\mathcal{F}}$
(3)

The free spectral range of the cavity is $\nu_{FSR}$, $\alpha_{chirp}$ is the
chirp rate, $\mathcal{F}$ is the finesse, and $L_{RT}$ is the round trip
optical path length of the cavity, twice its physical length for standing wave
geometries. At $T_{max}$ the Rabi frequency is reduced by a factor of two.
Beam splitter pulses at this Doppler shift must be twice as long as on
resonance and the velocity class that can be addressed is also reduced by a
factor of two. The sensitivity of circulating intensity to perturbations in
laser frequency or cavity length are all dramatically increased as the laser
moves off of resonance. The first order insensitivity on resonance is replaced
by the steep gradient of the Lorentzian Cavity transmission function (?).

### Pockels Cell

The Pockels cell used is a Gooch and Housego Light Gate 3 BBO Pockels cell, AR
coated at 780 nm. Precautions should be taken to avoid damage to the crystal
caused by the sustained application of high voltages (?).

### Laser system

Light at 780 nm, resonant with the D2 line of rubidium 87, is generated in
fiber by frequency doubling a 1560 nm telecoms laser as shown in figure 4A.
The output from the 1560 nm seed, a Koheras BASIK E15 from NKT, is split into
two arms for the two Raman frequencies. One of the arms is also used to
stabilize the laser to the cavity using the Pound-Drever-Hall (PDH) technique.
In both arms light travels through an electro-optic modulator (EOM), MPZ-LN-10
from iXblue, is amplified by an Erbium-doped fiber amplifier (EDFA), a CEFA-C-
PB from Keopsys, frequency doubled by a quasi-phase matched Lithium-Niobate
waveguide for second harmonic generation (SHG), WH-0780 from NTT Electronics,
and modulated using an acousto-optic modulator (AOM), MT110-NIR20 from AA
Opto-Electronic, before passing through an OZ Optics optical isolator to a
fiber collimator where the light is coupled into free space.

The EOM on the arm not locked to the cavity is modulated at 6.835 GHz; the
ground state hyperfine splitting in rubidium 87. To produce this a 7 GHz clock
is mixed with 164 MHz from an ARTIQ DDS channel using a single sideband mixer
achieving 32.7 dB carrier suppression. The RF signals for the remaining EOM
and AOMs are supplied by other ARTIQ DDS channels and an arbitrary function
generator.

### Locking

As we inject light into both cavity polarizations, the standard locking scheme
using a polarizing beamsplitter and a $\lambda/4$ (quarter-wave plate) to
separate the incident and reflected beams to be sampled for PDH stabilization
is not viable (?). Instead we insert a 5% beam sampler into one of the two
polarization beam paths, and stabilize the laser frequency on the sampled
fraction of the light reflected from the cavity with standard PDH. This loop
is run continuously at low optical power and a large detunig to maintain laser
lock whilst minimizing spontaneous emission. The laser is pulsed to full power
for the Raman pulses. Adding a sampler and photodiode to the second
polarization would allow the Pockels cell voltage to be actively stabilized
rather than the existing open-loop configuration.

### Cavity mode matching

Efficient mode matching is achieved by producing good spatial overlap between
the input beam and the cavity mode. We designed a custom beam delivery board
to overlap the two laser polarizations on a polarizing beamsplitter, before
passing through a telescope, $f_{3}$ and $f_{4}$ on figure 4 to produce a mode
of the correct size on the input mirror. Mechanical slides allow the telescope
lenses to be translated along the optical axis without introducing tilt or
transverse misalignments. The input diameter and location is measured with a
CCD which we are able to translate along the cavity rail around the location
of the input mirror (removed for this procedure). The telescope lenses are
adjusted until the target diameter, determined analytically, is achieved.

Figure 4: (A): Laser system schematic. Light at 780 nm is produced by
amplifying and frequency doubling a 1560 nm seed. (B): Cavity beam delivery
board. Mirrors overlap orthogonal linear polarized light from two fiber
collimators. The light is reflected through a mode matching telescope and
mirrors $m_{3}$ and $m_{4}$ are used to achieve spatial overlap with the
cavity mode. The cavity reflected light from the polarization shown in blue is
reflected by a 5 % beam sampler onto a bespoke photodiode for locking using
the Pound-Drever-Hall technique.

Spatial overlap with the cavity mode is attained with a final mirror pair on
the delivery board $m_{3}$ and $m_{4}$. These mirrors are adjusted to maximize
coupling to the fundamental mode (and minimize coupling to higher orders) as
observed on the cavity transmission photodiode (see figure 3). Precise
adjustment of these optics allows mode matching efficiencies above 90% before
adding the Pockels cell.

### Optics alignment

Despite the large number of optics present in this cavity, optical alignment
is relatively straightforward. Beginning with an optic axis parallel to the
rail, optics are added in turn and adjusted to overlap back reflections with
the input beams, see figure 1A. Coarse spacing is measured with mm-level
precision, and the critical telescope spacing is adjusted with a micrometer
stage to deliver the target beam size from the fixed input, as measured at the
waist of the large mode region. Retro reflection from $m_{4}$ in the absence
of $m_{3}$ allows the orientation of $m_{4}$ to be precisely defined. $m_{3}$
is then added, and coarsely aligned by making similar adjustments. Final
adjustments should be completed by optimizing the transmission signal. The
Pockels cell requires a well centered beam and accurate alignment of roll, tip
and tilt, carried out once the cavity is fully aligned. The coarse alignment
procedure follows that of standard q-switching applications (?), whilst the
final fine tuning is best achieved by supplying a triangle wave voltage to the
cell and adjusting the Pockels cell tip and tilt to minimize voltage induced
coupling to higher order modes as observed on the cavity transmission.
Temporarily adding a half waveplate on the input to the cavity simplifies
precise roll alignment, and enables projection of a single input polarization
on both crystal axis for alignment purposes. When the crystal axis is aligned
with the cavity axis, the applied voltage induces birefringence, splitting the
transmission peak in two and increasing this separation with increased
voltage, see figure 3. The appearance of other peaks or variations in peak
amplitude indicate misalignment.

## References

  * 1. K. Bongs, M. Holynski, J. Vovrosh, P. Bouyer, G. Condon, E. Rasel, C. Schubert, W. P. Schleich, A. Roura, Taking atom interferometric quantum sensors from the laboratory to real-world applications.  Nature Reviews Physics pp. 1–9 (2019).
  * 2. C. Freier, M. Hauth, V. Schkolnik, B. Leykauf, M. Schilling, H. Wziontek, H.-G. Scherneck, J. Müller, A. Peters, Mobile quantum gravity sensor with unprecedented stability.  Journal of Physics: Conference Series 723, 012050 (2016).
  * 3. B. Barrett, R. Geiger, I. Dutta, M. Meunier, B. Canuel, A. Gauguet, P. Bouyer, A. Landragin, The sagnac effect: 20 years of development in matter-wave interferometry.  Comptes Rendus Physique 15, 875 - 883 (2014). The Sagnac effect: 100 years later / L’effet Sagnac : 100 ans après.
  * 4. D. Budker, M. Romalis, Optical magnetometry.  Nature Physics 3, 227-234 (2007).
  * 5. S. Bize, The unit of time: Present and future directions.  Comptes Rendus Physique 20, 153 - 168 (2019). The new International System of Units / Le nouveau Système international d’unités.
  * 6. A. Trimeche, M. Langlois, S. Merlet, F. Pereira Dos Santos, Active control of laser wavefronts in atom interferometers.  Phys. Rev. Applied 7, 034016 (2017).
  * 7. R. Karcher, A. Imanaliev, S. Merlet, F. P. Dos Santos, Improving the accuracy of atom interferometers with ultracold sources.  New Journal of Physics 20, 113041 (2018).
  * 8. M.-k. Zhou, Q. Luo, L.-l. Chen, X.-c. Duan, Z.-k. Hu, Observing the effect of wave-front aberrations in an atom interferometer by modulating the diameter of raman beams.  Phys. Rev. A 93, 043610 (2016).
  * 9. V. Schkolnik, B. Leykauf, M. Hauth, C. Freier, A. Peters, The effect of wavefront aberrations in atom interferometry.  Applied Physics B 120, 311–316 (2015).
  * 10. J. M. Hogan, D. M. Johnson, S. Dickerson, T. Kovachy, A. Sugarbaker, S.-w. Chiow, P. W. Graham, M. A. Kasevich, B. Saif, S. Rajendran, et al., An atomic gravitational wave interferometric sensor in low earth orbit (agis-leo).  General Relativity and Gravitation 43, 1953–2009 (2011).
  * 11. L. Zhou, Z. Y. Xiong, W. Yang, B. Tang, W. C. Peng, K. Hao, R. B. Li, M. Liu, J. Wang, M. S. Zhan, Development of an atom gravimeter and status of the 10-meter atom interferometer for precision gravity measurement.  General Relativity and Gravitation 43, 1931-1942 (2011).
  * 12. S. wey Chiow, T. Kovachy, J. M. Hogan, M. A. Kasevich, Generation of 43 w of quasi-continuous 780 nm laser light via high-efficiency, single-pass frequency doubling in periodically poled lithium niobate crystals.  Opt. Lett. 37, 3861–3863 (2012).
  * 13. P. Hamilton, M. Jaffe, J. M. Brown, L. Maisenbacher, B. Estey, H. Müller, Atom interferometry in an optical cavity.  Phys. Rev. Lett. 114, 100405 (2015).
  * 14. I. Riou, N. Mielec, G. Lefèvre, M. Prevedelli, A. Landragin, P. Bouyer, A. Bertoldi, R. Geiger, B. Canuel, A marginally stable optical resonator for enhanced atom interferometry.  Journal of Physics B: Atomic, Molecular and Optical Physics 50, 155002 (2017).
  * 15. P. W. Graham, J. M. Hogan, M. A. Kasevich, S. Rajendran, R. W. Romani, Mid-band gravitational wave detection with precision atomic sensors.  arXiv preprint arXiv:1711.02225 (2017).
  * 16. L. Badurina, E. Bentine, D. Blas, K. Bongs, D. Bortoletto, T. Bowcock, K. Bridges, W. Bowden, O. Buchmueller, C. Burrage, J. Coleman, G. Elertas, J. Ellis, C. Foot, V. Gibson, M. G. Haehnelt, T. Harte, S. Hedges, R. Hobson, M. Holynski, T. Jones, M. Langlois, S. Lellouch, M. Lewicki, R. Maiolino, P. Majewski, S. Malik, J. March-Russell, C. McCabe, D. Newbold, B. Sauer, U. Schneider, I. Shipsey, Y. Singh, M. A. Uchida, T. Valenzuela, M. van der Grinten, V. Vaskonen, J. Vossebeld, D. Weatherill, I. Wilmut, Aion: An atom interferometer observatory and network (2019).
  * 17. D. O. Sabulsky, I. Dutta, E. A. Hinds, B. Elder, C. Burrage, E. J. Copeland, Experiment to detect dark energy forces using atom interferometry.  Phys. Rev. Lett. 123, 061102 (2019).
  * 18. G. Lamporesi, A. Bertoldi, L. Cacciapuoti, M. Prevedelli, G. M. Tino, Determination of the newtonian gravitational constant using atom interferometry.  Phys. Rev. Lett. 100, 050801 (2008).
  * 19. R. H. Parker, C. Yu, W. Zhong, B. Estey, H. Müller, Measurement of the fine-structure constant as a test of the standard model.  Science 360, 191–195 (2018).
  * 20. C. Overstreet, P. Asenbaum, T. Kovachy, R. Notermans, J. M. Hogan, M. A. Kasevich, Effective inertial frame in an atom interferometric test of the equivalence principle.  Phys. Rev. Lett. 120, 183604 (2018).
  * 21. P. Asenbaum, C. Overstreet, M. Kim, J. Curti, M. A. Kasevich, Atom-interferometric test of the equivalence principle at the ${10}^{-12}$ level.  Phys. Rev. Lett. 125, 191101 (2020).
  * 22. S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kasevich, S. Rajendran, Atomic gravitational wave interferometric sensor.  Phys. Rev. D 78, 122002 (2008).
  * 23. M. Jaffe, P. Haslinger, V. Xu, P. Hamilton, A. Upadhye, B. Elder, J. Khoury, H. Müller, Testing sub-gravitational forces on atoms from a miniature in-vacuum source mass.  Nature Physics 13, 938–942 (2017).
  * 24. X. Wu, Z. Pagel, B. S. Malek, T. H. Nguyen, F. Zi, D. S. Scheirer, H. Müller, Gravity surveys using a mobile atom interferometer.  Science advances 5, eaax0800 (2019).
  * 25. A. Hinton, M. Perea-Ortiz, J. Winch, J. Briggs, S. Freer, D. Moustoukas, S. Powell-Gill, C. Squire, A. Lamb, C. Rammeloo, B. Stray, G. Voulazeris, L. Zhu, A. Kaushik, Y.-H. Lien, A. Niggebaum, A. Rodgers, A. Stabrawa, D. Boddice, S. R. Plant, G. W. Tuckwell, K. Bongs, N. Metje, M. Holynski, A portable magneto-optical trap with prospects for atom interferometry in civil engineering.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20160238 (2017).
  * 26. Y. Bidel, N. Zahzam, A. Bresson, C. Blanchard, M. Cadoret, A. V. Olesen, R. Forsberg, Absolute airborne gravimetry with a cold atom sensor.  Journal of Geodesy 94, 20 (2020).
  * 27. M. Dovale-Álvarez, D. D. Brown, A. W. Jones, C. M. Mow-Lowry, H. Miao, A. Freise, Fundamental limitations of cavity-assisted atom interferometry.  Phys. Rev. A 96, 053820 (2017).
  * 28. M. D. Álvarez, Optical cavities for optical atomic clocks, atom interferometry and gravitational-wave detection (Springer, 2019).
  * 29. S. L. Kristensen, M. Jaffe, V. Xu, C. D. Panda, H. Müller, Raman transitions driven by phase-modulated light in a cavity atom interferometer (2020).
  * 30. H. Kogelnik, T. Li, Laser beams and resonators.  Appl. Opt. 5, 1550–1567 (1966).
  * 31. A. Siegman, Lasers (University Science Books, 1986).
  * 32. S. M. Dickerson, J. M. Hogan, A. Sugarbaker, D. M. S. Johnson, M. A. Kasevich, Multiaxis inertial sensing with long-time point source atom interferometry.  Phys. Rev. Lett. 111, 083001 (2013).
  * 33. G. Ohio, Pockels cell alignment for single pass systems (2012).
  * 34. E. D. Black, An introduction to pound–drever–hall laser frequency stabilization.  American Journal of Physics 69, 79-87 (2001).

Acknowledgements: The authors thank Camron Nourshargh, Farzad Hayati, and Luuk
Earl for their contributions.
Funding: This work was supported by the EPSRC under grant number EP/T001046/1
as part of the UK National Quantum Technologies Programme.
Competing Interests: The authors declare that they have no competing
interests.
Data and materials availability: Additional data related to this paper may be
requested from the authors.


EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 LHCb-DP-2021-001 May 17, 2021

A parametrized Kalman filter for fast track fitting at LHCb

P. Billoir1, M. De Cian2, P. A. Günther3, S. Stemmle3,†

1LPNHE, Sorbonne Université, Paris Diderot Sorbonne Paris Cité, CNRS/IN2P3,
Paris, France 2Institute of Physics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), Lausanne, Switzerland
3Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
†Author was at institute at time work was performed.

We present an alternative implementation of the Kalman filter employed for
track fitting within the LHCb experiment. It uses simple parametrizations for
the extrapolation of particle trajectories in the field of the LHCb dipole
magnet and for the effects of multiple scattering in the detector material. A
speedup of more than a factor of four is achieved while maintaining the
quality of the estimated track quantities. This Kalman filter implementation
could be used in the purely software-based trigger of the LHCb upgrade.

Published in Computer Physics Communications 265, 108026 (2021)

© 2024 CERN for the benefit of the LHCb collaboration. CC-BY-4.0 licence.

##  1 Introduction

The LHCb experiment is a dedicated heavy flavour physics experiment at the LHC
focusing on the study of hadrons containing $b$ and $c$ quarks [1]. Due to the
high luminosity at the LHC and the high proton-proton interaction cross
section, a sophisticated trigger system is needed to reduce the rate of
collisions saved for offline analysis. During Runs 1 and 2 of the LHC, this
trigger system consisted of a hardware stage, reducing the rate from
$40\text{\,MHz}$ to $1\text{\,MHz}$, followed by a two-stage software trigger.
In the latter, the full tracking system was read out and a partial (first
stage) and full (second stage) event reconstruction were performed [2]. Both
software stages included a fit of selected track candidates using a Kalman
filter to extract their parameters and to reject fake tracks. In addition, the
software trigger allowed an online calibration and alignment of the detector
[3].

During Run 3 of the LHC, LHCb will be provided with a factor five higher
luminosity compared to Run 2. In this scope, most of the subdetectors are
currently being replaced or upgraded[4, 5, 6, 7] and a new trigger strategy
has been developed [8]. The hardware trigger will be removed and a two-stage,
fully software-based trigger will process the full $30\text{\,MHz}$111The
nominal bunch-crossing frequency of the LHC is 40 MHz, however empty and non-
colliding bunches reduce this to a collision frequency of 30 MHz at LHCb. of
bunch-crossing rate. In the first stage, tracks with a high transverse
momentum ($p_{\mathrm{T}}$) and primary vertices will be reconstructed. These
objects are used to select events with displaced topologies typical for
$b$-hadron and $c$-hadron decays, and to select high-$p_{\mathrm{T}}$ objects
from decays of heavy vector bosons. In the second stage, a full event
reconstruction will be performed, without any requirement on the
$p_{\mathrm{T}}$ and including particle identification. A large number of
exclusive and several universal event selections based on the decay topology
will be applied.

In LHCb, track reconstruction is split into a pattern recognition and a Kalman
filtering[9, 10] stage. During pattern recognition, sets in each subdetector
are constructed from signals that potentially result from the passage of a
single charged particle. Simple parametrizations are used throughout this
procedure as it is only concerned with finding the right sets of signals and
not to provide the best estimate of the track parameters. During the filtering
stage, an estimate for the track parameters is calculated, and fake tracks are
rejected. Given that the output of the filtering stage is used for physics
selections the best possible precision needs to be achieved, hence an
(extended) Kalman filter is used for track fitting. Ideally, Kalman filtering
of the track candidates is already performed during the first trigger stage.
However, the Kalman filter which was used during Run 1 and 2 in LHCb, in the
following called default Kalman, is significantly too slow. It relies on
lookup tables for the magnetic field and the material distribution of the
detector[11], so-called maps. In addition it uses Runge-Kutta methods to solve
the differential equations necessary to propagate the particle through the
regions with an inhomogeneous magnetic field. Accessing the values in the
lookup table and solving the differential equations are time consuming and
prohibit the usage of the current Kalman filter in the first stage of the
upgraded trigger system. This conclusion is independent of the choice of
computing architecture (CPU or GPU) which is used for the first trigger stage.

In this paper, a fully parametrized version of the Kalman filter in LHCb,
called parametrized Kalman, is presented. It obtains precise values of track
parameters and track quality variables, while relying on neither
computationally costly extrapolation methods nor material or magnetic field
maps.

##  2 Detector and simulation

The LHCb detector [1] is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$. Its Run 3 configuration includes a high-
precision tracking system consisting of a silicon-pixel vertex detector
surrounding the $pp$ interaction region [5] (VELO), a large-area silicon-strip
detector (Upstream Tracker (UT)) [7] located upstream of a dipole magnet with
a bending power of about $4{\mathrm{\,Tm}}$ [12], and three stations of
scintillating-fibre detectors (SciFi) [7] placed downstream of the magnet.
Different types of charged hadrons are distinguished using information from
two ring-imaging Cherenkov detectors [13, 6]. Photons, electrons and hadrons
are identified by a calorimeter system consisting of an electromagnetic and a
hadronic calorimeter [14, 6]. Muons are identified by a system composed of
alternating layers of iron and multiwire proportional chambers [15, 6].

Given the lack of collision data at this point for Run 3, simulation is
required to model the effects of the detector response, the detector
acceptance and the imposed selection requirements. In the simulation, $pp$
collisions are generated using Pythia [16, *Sjostrand:2006za] with a specific
LHCb configuration [18]. Decays of unstable particles are described by EvtGen
[19], in which final-state radiation is generated using Photos [20]. The
interaction of the generated particles with the detector, and its response,
are implemented using the Geant4 toolkit [21, *Agostinelli:2002hh] as
described in Ref. [23].

##  3 Principles

In the following, the Kalman filter formalism and its application in the LHCb
track reconstruction is outlined. During Kalman filtering, the information
from measurements at detector planes is successively combined to obtain
optimal estimates of the track parameters. The track is represented as a set
of states at fixed $z$-positions222The detector coordinate system is chosen
such that the $z$-axis is parallel to the beam line and charged particles are
deflected in the direction of the $x$-axis., which are typically detector
layers. Each of these states is given by
$\boldsymbol{x}=(x,y,t_{x},t_{y},\frac{q}{p})$ and the corresponding
covariance matrix $\boldsymbol{P}$, where $t_{x}$ and $t_{y}$ are the slopes
with respect to the $z$ axis, $q$ the charge of the particle in units of the
electron charge and $p$ its absolute momentum.

The Kalman filter procedure needs an estimate of a state as a starting point.
Filtering is then a repeated application of two steps. Firstly, the current
state is extrapolated to the next detector layer, and secondly, the
extrapolated state is updated using the measurement in this layer. If the
track has no associated measurement in this layer, the update step is omitted.
These steps can be formalized as follows: given the state
($\boldsymbol{x}_{k-1|k-1}$, $\boldsymbol{P}_{k-1|k-1}$) at position
$z_{k-1}$, the extrapolated state ($\boldsymbol{x}_{k|k-1}$,
$\boldsymbol{P}_{k|k-1}$) at position $z_{k}$ is given by

$\displaystyle\boldsymbol{x}_{k|k-1}$
$\displaystyle=\boldsymbol{f}_{k}(\boldsymbol{x}_{k-1|k-1}),$ (1)
$\displaystyle\boldsymbol{P}_{k|k-1}$
$\displaystyle=\boldsymbol{F}_{k}\boldsymbol{P}_{k-1|k-1}\boldsymbol{F}_{k}^{T}+\boldsymbol{Q}_{k},$
(2)

where the extrapolation function $\boldsymbol{f}_{k}(\boldsymbol{x})$ is given
by five individual mappings
$\boldsymbol{f}_{k}=(f_{k}^{x},f_{k}^{y},f_{k}^{t_{x}},f_{k}^{t_{y}},f_{k}^{\frac{q}{p}})$.
This leads to the transport matrix $\boldsymbol{F}_{k}$ as

$\displaystyle F_{k}^{ij}=\frac{\partial f_{k}^{i}}{\partial x_{j}}.$ (3)

The noise matrix $\boldsymbol{Q}_{k}$ accounts for uncertainties of the
extrapolation, e.g. due to scattering at the material of the detector layers
or the material in between.

The extrapolated state is then combined with the measurement
$\boldsymbol{m}_{k}$ in the respective detector layer to obtain the new state
estimate at the position $\boldsymbol{z}_{k}$, $\boldsymbol{x}_{k|k}$ and
$\boldsymbol{P}_{k|k}$, using the following steps:

$\displaystyle\boldsymbol{r}_{k}$
$\displaystyle=\boldsymbol{m}_{k}-\boldsymbol{H}_{k}\boldsymbol{x}_{k|k-1},$
(4) $\displaystyle\boldsymbol{S}_{k}$
$\displaystyle=\boldsymbol{H}_{k}\boldsymbol{P}_{k|k-1}\boldsymbol{H}_{k}^{T}+\boldsymbol{R}_{k},$
(5) $\displaystyle\boldsymbol{K}_{k}$
$\displaystyle=\boldsymbol{P}_{k|k-1}\boldsymbol{H}_{k}^{T}\boldsymbol{S}_{k}^{-1},$
(6) $\displaystyle\boldsymbol{x}_{k|k}$
$\displaystyle=\boldsymbol{x}_{k|k-1}+\boldsymbol{K}_{k}\boldsymbol{r}_{k},$
(7) $\displaystyle\boldsymbol{P}_{k|k}$
$\displaystyle=(\boldsymbol{1}-\boldsymbol{K}_{k}\boldsymbol{H}_{k})\boldsymbol{P}_{k|k-1}.$
(8)

Here $\boldsymbol{H}_{k}$ projects the estimated state vector to the
measurement space in order to allow a calculation of the residual
$\boldsymbol{r}_{k}$. The covariance matrix of this residual is given by
$\boldsymbol{S}_{k}$ and is combined with the covariance matrix of the state
to obtain the Kalman gain $\boldsymbol{K}_{k}$. The latter defines then how
the estimated state is modified by the residual. The variance of the residual
is given by $\boldsymbol{R}_{k}$.

Starting at the most upstream measurement, the measurements are successively
added and the track parameters updated until the last detector layer is
reached. The same procedure is repeated starting at the most downstream
measurement and successively including more upstream measurements. This yields
two sets of states at every measurement position, which can be combined to
obtain the respective optimal state.

The quality of a track can be estimated by its $\chi^{2}_{\text{track}}$
value. The value at each measurement is given by:

$\displaystyle\chi^{2}_{k}$
$\displaystyle=\chi^{2}_{k-1}+\boldsymbol{r}_{k}^{T}\boldsymbol{P}_{k|k}^{-1}\boldsymbol{r}_{k},$
(9)

and $\chi^{2}_{\text{track}}$ is then simply $\chi^{2}_{k}$ after all
measurements have been added using the combined, optimal states.

The optimal state estimates and the measurement information can also be used
to remove measurements that show a large separation from the fitted trajectory
by having a large contribution to the $\chi^{2}_{\text{track}}$ value. They
are therefore likely to be wrongly associated to the respective track, and are
so-called outliers. Once an outlier is removed, all Kalman filter steps are
performed again. This procedure can be repeated until the maximum allowed
number of outliers are removed, or no more outliers are present.

The above formalism is also the basis of the Kalman filter that is currently
used for track fitting in the LHCb experiment. The extrapolation functions
$\boldsymbol{f}_{k}$ are based on maps of the magnetic field along the
trajectory and numerical models for the extrapolations. Their complexities
range up to a fifth-order Runge-Kutta method. The noise matrices
$\boldsymbol{Q_{k}}$ are obtained by a dedicated model for the multiple
scattering and a map of the material traversed by the particle.

In the parametrized Kalman filter presented in this paper, these two costly
steps are replaced by simple parametrizations. The extrapolation functions
$\boldsymbol{f}_{k}$ are given by analytic expressions that allow a fast
evaluation and calculation of the derivatives in Equation 3. The noise
matrices $\boldsymbol{Q}_{k}$ depend on the momentum of the particle and are
parametrized by a few parameters per extrapolation step.

An important difference with respect to the default Kalman filter is the
treatment of energy loss due to the interaction with the detector material.
While the multiple scattering is taken directly into account, the energy loss
is not part of the extrapolation functions $\boldsymbol{f}_{k}$, i.e.
$f_{k}^{\frac{q}{p}}$ is the unity transformation. This shortcoming is
compensated by choosing the momentum of the state vectors to represent the
momentum at the moment of production of the particle. Thereby, the
extrapolation functions also take this initial momentum as input and thus
indirectly take into account all energy loss that happened on average up to
the respective detector layer. The only caveat being that $\frac{q}{p}$ after
the filtering is only the best representation of the true value at the
production point of the particle.

##  4 Parametrizations

Depending on the strength of the magnetic field and the typical distance
between detector layers, different empirical analytical functions for the
extrapolation are used.

Inside the VELO, where the magnetic field is very weak, these functions and
the noise matrix are given by:

$\displaystyle\boldsymbol{f}(\boldsymbol{x})=\begin{pmatrix}f^{x}(\boldsymbol{x})\\\
f^{y}(\boldsymbol{x})\\\ f^{t_{x}}(\boldsymbol{x})\\\
f^{t_{y}}(\boldsymbol{x})\\\
f^{\frac{q}{p}}(\boldsymbol{x})\end{pmatrix}=\begin{pmatrix}x+0.5[t_{x}+f^{t_{x}}(\boldsymbol{x})]\Delta
z\ \\\ y+t_{y}\Delta z\\\
t_{x}+p^{\text{V}}_{0}\frac{q}{p}(z_{0}+p^{\text{V}}_{1})\Delta z\\\ t_{y}\\\
\frac{q}{p}\end{pmatrix}$ (10)

and

$\displaystyle\boldsymbol{Q}=\begin{pmatrix}\left(\tilde{p}^{\text{V}}_{1}\Delta
z\right)^{2}Q^{t_{x}t_{x}}&0&\tilde{p}^{\text{V}}_{2}\sqrt{Q^{xx}Q^{t_{x}t_{x}}}&0&0\\\
0&\left(\tilde{p}^{\text{V}}_{1}\Delta
z\right)^{2}Q^{t_{y}t_{y}}&0&\tilde{p}^{\text{V}}_{3}\sqrt{Q^{yy}Q^{t_{y}t_{y}}}&0\\\
\tilde{p}^{\text{V}}_{2}\sqrt{Q^{xx}Q^{t_{x}t_{x}}}&0&\left(\tilde{p}^{\text{V}}_{0}\left|\frac{q}{p}\right|\right)^{2}&0&0\\\
0&\tilde{p}^{\text{V}}_{3}\sqrt{Q^{yy}Q^{t_{y}t_{y}}}&0&\left(\tilde{p}^{\text{V}}_{0}\left|\frac{q}{p}\right|\right)^{2}&0\\\
0&0&0&0&0\\\ \end{pmatrix},$ (11)

where $\Delta z$ is the extrapolation distance along the $z$-direction and
$z_{0}$ the initial or final $z$ coordinate for a downstream or upstream
extrapolation, respectively. The parameters $p^{\text{V}}_{0}$,
$p^{\text{V}}_{1}$ and $\tilde{p}^{\text{V}}_{0}$ to
$\tilde{p}^{\text{V}}_{3}$ are the same for all upstream and downstream
extrapolations inside the VELO. They are determined using simulated
${{B}^{0}_{s}}\\!\rightarrow\phi\phi$ decays within the LHCb software
framework, where $\phi\\!\rightarrow{{K}^{+}}{{K}^{-}}$. This simulated sample
allows to create a dataset $D$, containing pairs of states representing two
consecutive measurements of one track inside the VELO. In addition to the true
state parameters obtained from the simulation, also an extrapolation of each
state to the $z$ position of the respective other state is included in the
dataset. Such extrapolation is based on the default extrapolation algorithm in
LHCb [11]. This dataset allows tuning the parameters employing a minimization
of the following likelihood-inspired function:

$\displaystyle\prod_{D}\left[\mathcal{G}\left(f^{s}(\boldsymbol{x_{1}})-\boldsymbol{x_{2}}^{s},\sqrt{Q^{ss}}\right)+c\right].$
(12)

Here, $\mathcal{G}(x,\sigma_{x})$ is a normalized Gaussian distribution
centered around $0$ with width $\sigma_{x}$. The two states of each dataset
entry are represented by $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$, and
the variable $s$ is one of the state variables, $s\in\\{x,t_{x},y,t_{y}\\}$.
The positive empirical constant $c$ is chosen to be small with respect to the
amplitude of the Gaussian function and softens the impact of outliers.

In a first step, the extrapolation functions $f^{x}$ to $f^{t_{x}}$ are tuned
individually, taking into account that $f^{x}$ depends on the previously
determined parameters for $f^{t_{x}}$. These tuning minimizations employ the
state vector $\boldsymbol{x}_{2}$ that is obtained by the extrapolation of the
state vector $\boldsymbol{x}_{1}$. This choice improves the precision of the
parametrized extrapolation, by removing the effect of multiple scattering that
would be present if instead the true state was chosen for
$\boldsymbol{x}_{2}$.

In a second step, the parameters of the extrapolation functions are fixed, and
a minimization of the following function is performed:

$\displaystyle\prod_{D}\left[\mathcal{G}_{2}\left(f^{d}(\boldsymbol{x_{1}})-\boldsymbol{x_{2}}^{d},f^{t_{d}}(\boldsymbol{x_{1}})-\boldsymbol{x_{2}}^{t_{d}},\sqrt{Q^{dd}},\sqrt{Q^{t_{d}t_{d}}},Q^{dt_{d}}/\sqrt{Q^{dd}Q^{t_{d}t_{d}}}\right)+c\right].$
(13)

Here, $\mathcal{G}_{2}(x,y,\sigma_{y},\sigma_{y},\rho)$ is a normalized two-
dimensional Gaussian distribution centered around $0$ with widths $\sigma_{x}$
and $\sigma_{y}$ and a correlation factor $\rho$. The variable $d$ is either
$x$ or $y$. In this minimization, the true state vector $\boldsymbol{x}_{2}$
is used in order to get the correct estimate of the parameters for the
respective elements of the noise matrix $\boldsymbol{Q}$.

Inside the UT and the SciFi detector stations, the magnetic field is
significantly stronger than inside the VELO and higher order terms are needed
for the extrapolation functions:

$\displaystyle\boldsymbol{f}(\boldsymbol{x})=\begin{pmatrix}x+\left[p^{\text{T}}_{3}t_{x}+(1-p^{\text{T}}_{3})f^{t_{x}}(\boldsymbol{x})\right]\Delta
z\\\
y+\left[p^{\text{T}}_{5}t_{y}+(1-p^{\text{T}}_{5})f^{t_{y}}(\boldsymbol{x})\right]\Delta
z\\\
t_{x}+\left[p^{\text{T}}_{0}\frac{q}{p}+p^{\text{T}}_{1}(\frac{q}{p})^{3}+p^{\text{T}}_{2}y^{2}\frac{q}{p}\right]\Delta
z\\\ t_{y}+p^{\text{T}}_{4}\frac{q}{p}t_{x}\frac{y}{|y|}\\\
\frac{q}{p}\end{pmatrix}.$ (14)

The noise matrix is given in full analogy to Equation 11 with the parameters
$\tilde{p}^{\text{T}}_{0}$ to $\tilde{p}^{\text{T}}_{3}$, where T either
stands for the UT or the SciFi detector. These parameters and the parameters
$p^{\text{T}}_{0}$ to $p^{\text{T}}_{4}$ are individually determined on
simulation for every step from one detector layer to the next and for the
upstream and downstream extrapolation separately. The same strategy as for the
tuning of the parameters related to the extrapolation inside the VELO is
followed.

For the long extrapolations between the different tracking subdetectors, more
sophisticated parametrizations are necessary. In the case of the step between
the VELO and the UT, where the magnetic field is still weak, the extrapolation
is based on two equations. The first describes the change in momentum along
the $x$-direction of the particle:

$\displaystyle\Delta
p_{x}=p\left(\frac{t_{x,\text{UT}}}{\sqrt{1+t^{2}_{x,\text{UT}}+t^{2}_{y,\text{UT}}}}-\frac{t_{x,\text{V}}}{\sqrt{1+t^{2}_{x,\text{V}}+t^{2}_{y,\text{V}}}}\right)=q\int\left(\text{d}\boldsymbol{l}\times\boldsymbol{B}\right)_{x},$
(15)

where $t_{x/y,\text{UT}}$ and $t_{x/y,\text{V}}$ are the state variables at
the first UT detector layer and the last measurement inside the VELO,
respectively. The right hand side of the equation consists of an integral of
the magnetic field along the trajectory of the particle. Note that the
integral expression is simply a parameter which was fitted for on the dataset.
The second ingredient for the extrapolation is to model the effect of the
magnetic field as a single kink of the trajectory at a certain $z$-position
$z_{\text{mag}}$ between the VELO and the UT:

$\displaystyle
x_{\text{UT}}=x_{\text{V}}+(z_{\text{mag}}-z_{\text{V}})t_{x,\text{V}}+(z_{\text{UT}}-z_{\text{mag}})t_{x,\text{UT}},$
(16)

where $z_{\text{V}}$ and $z_{\text{UT}}$ are the positions of the states
inside the VELO and the UT, respectively.

Equation 15 can be solved for $t_{x,\text{UT}}$ and Equation 16 is then
employed to get an expression for $x_{\text{UT}}$. The unknowns in these
expressions are parametrized as a function of the state variables inside the
VELO:

$\displaystyle t_{y,\text{UT}}$
$\displaystyle=t_{y,\text{V}}+p^{\text{S}}_{0}\frac{q}{p}t_{x,\text{V}}\frac{y_{\text{V}}}{|y_{\text{V}}|}$
(17)
$\displaystyle\int\left(\text{d}\boldsymbol{l}\times\boldsymbol{B}\right)_{x}$
$\displaystyle=p^{\text{S}}_{1}+p^{\text{S}}_{2}z_{\text{V}}+p^{\text{S}}_{3}t^{2}_{y,\text{V}}$
(18) $\displaystyle z_{\text{mag}}$
$\displaystyle=p^{\text{S}}_{4}+p^{\text{S}}_{5}z_{\text{V}}+p^{\text{S}}_{6}z^{2}_{\text{V}}+p^{\text{S}}_{7}t^{2}_{y,\text{V}}.$
(19)

In addition, the $y$-position of the extrapolated state is given by:

$\displaystyle
y_{\text{UT}}=y_{\text{V}}+\left[p^{\text{S}}_{8}t_{y,\text{V}}+(1-p^{\text{S}}_{8})t_{y,\text{UT}}\right]\Delta
z,$ (20)

where $\Delta z$ is defined as the difference between $z_{\text{UT}}$ and
$z_{\text{V}}$. The noise matrix is defined in analogy to Equation 11 with the
parameters $\tilde{p}^{\text{S}}_{0}$ to $\tilde{p}^{\text{S}}_{3}$. These
parameters and the parameters $p^{\text{S}}_{0}$ to $p^{\text{S}}_{8}$ are
individually determined for the upstream and downstream extrapolation. The
same strategy as for the tuning of the parameters related to the extrapolation
inside the VELO is followed.

The extrapolation from the UT to the SciFi detector is more delicate because
it is done over a distance of more than 5 meters through a strong magnetic
field. Moreover, this field is far from uniform - in particular, it varies
rapidly in the upper and lower regions, close to the magnet yoke. To ensure a
good quality of the global track fit, the error on the extrapolation should be
well below the other sources of error, mainly multiple scattering. The chosen
solution is an expansion of the magnetic deviation in powers of $q/p$. The
parametrization aims at giving good precision for charged particles used in
physics analyses, that is for trajectories which roughly come from the origin.

To do so, the ideal direction $(t_{x}^{0},t_{y}^{0})$ as the one of a particle
of charge $q$, momentum $p$, starting from the origin and hitting the UT
detector layer in a given point $(x,y)$ is defined. As a good approximation,
we can take $t_{x}^{0}=x/z+{\cal B}q/p$, $t_{y}^{0}=y/z$, where $\cal{B}$ is
proportional to the integrated field between the origin and the UT. The
deviations from the ideal direction, $\delta t_{x}=t_{x}-t_{x}^{0}$, $\delta
t_{y}=t_{y}-t_{y}^{0}$, are small, so only a first order expansion in $\delta
t_{x},\delta t_{y}$ is considered. Corrections of higher order would be
negligible compared to multiple scattering errors.

Finally, a polynomial expansion in $q/p$ for the ideal direction is built, and
a correction in $\delta t_{x},\delta t_{y}$ with coefficients which are
themselves polynomials of $q/p$ is added:

$\displaystyle f^{x}(\boldsymbol{x})=x+t_{x}\Delta
z+\sum_{k=1}^{K_{1}}A^{x}_{k}(x,y)\left(\frac{q}{p}\right)^{k}+\sum_{k=1}^{K_{2}}\left(B^{x}_{k}(x,y)\,\delta
t_{x}+C^{x}_{k}(x,y)\,\delta t_{y}\right)\left(\frac{q}{p}\right)^{k},$ (21)

where the first two terms are the straight line extrapolation, and the next
ones the curvature correction. Similar expressions are used for the other
state parameters $f^{y}(\boldsymbol{x})$, $f^{t_{x}}(\boldsymbol{x})$,
$f^{t_{y}}(\boldsymbol{x})$. The degrees of expansion $K_{1}$ and $K_{2}$ are
tuned for each parameter to obtain the required precision. In practice
$K_{1}=9$, $K_{2}=7$ for $f^{x}$ and $f^{t_{x}}$ and $K_{1}=7$, $K_{2}=5$ for
$f^{y}$ and $f^{t_{y}}$ are used.

The dependence on $x,y$ of the coefficients $A^{u}_{k}$, $B^{u}_{k}$,
$C^{u}_{k}$, with $u=x,y,t_{x},t_{y}$, is described through a tabulation on a
grid of 50$\times$50 points regularly spaced on the rectangle defined by
$|x/z|\leq 0.25$, $|y/z|\leq 0.25$, by steps $\Delta X$, $\Delta Y$. In order
to avoid a systematic convexity bias of a bilinear interpolation, the values
at $x,y$ are computed by a quadratic interpolation between the tabulated
values at the six closest points on the grid: if $(X,Y)$ is the closest one,
these values are: $F_{00}=(X,Y)$, $F_{+0}=F(X+\Delta X,Y)$, $F_{-0}=F(X-\Delta
X,Y)$, $F_{0+}=F(X,Y+\Delta Y)$, $F_{0-}=F(X,Y-\Delta Y)$, and
$F_{\varepsilon_{x}\varepsilon_{y}}=F(X+\varepsilon_{x}\Delta
X,Y+\varepsilon_{y}\Delta Y)$, where $\varepsilon_{x}$ and $\varepsilon_{y}$
are the signs of $\xi=(x-X)/\Delta X$ and $\psi=(y-Y)/\Delta Y$, respectively.
With these notations the interpolation formula for a quantity $F$ is given by
:

$\displaystyle F(x,y)=F_{00}+F_{d}\,\xi\psi+\big{(}$
$\displaystyle(F_{+0}-F_{-0})\,\xi+(F_{0+}-F_{0-})\,\psi$
$\displaystyle+(F_{+0}+F_{-0}-2F_{00})\,\xi^{2}+(F_{0+}+F_{0-}-2F_{00})\,\psi^{2}\big{)}/2$
(22) $\displaystyle\text{with}\;\;\;F_{d}=$
$\displaystyle\varepsilon_{x}\varepsilon_{y}(F_{00}+F_{\varepsilon_{x}\varepsilon_{y}}-F_{\varepsilon_{x}0}-F_{0\varepsilon_{y}}).$
(23)

The tabulated values are obtained using the standard Runge-Kutta method of
order 4, with 20 values of $q/p$ in the range ($-1/p_{min},1/p_{min}$), with
$p_{min}=3000\text{\,Me\kern-1.00006ptV\\!/}c$ and a polynomial fit in $q/p$.
As a consequence, they do not give a reliable result for momenta below
$p_{min}$. Another limitation is the larger errors on the edges of the
acceptance, especially for $|t_{y}|\simeq 0.25$, where the field has strong
spatial variations.

##  5 Performance

A sample of simulated proton-proton collisions that include a
$B^{0}_{s}\rightarrow\phi\phi$, $\phi\\!\rightarrow{{K}^{+}}{{K}^{-}}$ decay
is used to compare the reconstruction quality of the parametrized and the
default Kalman filter. The extrapolation of the most upstream state estimate
to the beam line is the same in both filters and is based on a simplified
material map of the detector [11]. Therefore, not the state near the beam
line, but the state at the most upstream measurement is employed for the
comparison of the two Kalman filters. Although only tracks with measurements
in each of the subdetectors are considered for this study, this is in
principle not a requirement for operating the parameterized Kalman filter

Figure 1 compares the resolution of the momentum, the $x$-position and the
slope $t_{x}$ as a function of the true momentum of a particle.

Figure 1: Comparison of the resolution in simulation in (top left) momentum,
(top right) $x$-position and (bottom) slope $t_{x}$ between the default and
parametrized Kalman filter. The resolution is represented by the root mean
square of the residual distribution when comparing to the true value.

Since the position and slope are nearly exclusively determined by the
measurements in the VELO, where only a very weak magnetic field is present,
the parametrizations of the parametrized Kalman filter are sufficient to
obtain results comparable to the default Kalman filter in these variables. In
contrast, the momentum estimate strongly depends on the extrapolations in
regions with strong magnetic field. There, especially at momenta below
10$\text{\,Ge\kern-1.00006ptV\\!/}c$, an up to $20\%$ worse resolution is
observed for the parametrized Kalman filter.

The Kalman filter does not only provide an estimate of the state parameters,
but also a corresponding covariance matrix. In Figure 2 the pull distributions
of the estimated momentum, $x$-position and slope $t_{x}$ for the parametrized
Kalman filter are shown.

Figure 2: Pull distributions of the momentum, $x$-position and slope $t_{x}$
estimates of the parametrized Kalman filter at the most upstream measurement.
The given values correspond to the mean, width and root mean square of a
Gaussian function that is fitted to the distribution.

In all three cases, good uncertainty estimates are visible. However, in
analogy to the observations made for the resolution, the pull distribution of
the momentum features slightly more pronounced tails.

Besides the estimate of the state near the beam line, which is used for the
reconstruction of charged particles, an important output of the Kalman filter
is the fit quality described by the $\chi^{2}_{\text{track}}$ per degrees of
freedom $N_{\text{dof}}$. In Figure 3, this quantity is shown for the
parametrized Kalman filter for real tracks coming from a particle and fake
tracks consisting of random combinations of clusters. In addition, the real
track efficiencies and fake track rejection rates are shown for both Kalman
filter versions when applying upper bounds on this quantity.

Figure 3: Track quality estimate, $\chi^{2}_{\text{track}}/N_{\text{dof}}$, in
simulation for the parametrized filter (left). Fake tracks are shown in red
and real tracks in black. Real track efficiency and fake track rejection for
the parametrized and default Kalman filter (right).

The parametrized Kalman filter shows a slightly worse but overall comparable
performance in separating the two track classes.

The fitted tracks are combined to reconstruct $B^{0}_{s}\rightarrow\phi\phi$
candidates. Figure 4 shows the invariant mass distribution of candidates based
on the two Kalman filter versions.

Figure 4: Reconstructed $B^{0}_{s}$ mass in simulated
$B^{0}_{s}\rightarrow\phi\phi$ decays for the parametrized and the default
Kalman filter. Fit projections are overlaid.

A single Gaussian distribution and a first order polynomial are employed to
model the signal peak and the combinatorial background, respectively. This
yields nearly identical estimated mass resolutions of
$12.8\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ and
$12.9\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ for the default and the
parametrized Kalman filter, respectively.

In order to compare the timing performance of the parametrized Kalman filter
and the default Kalman filter, throughput studies on a machine with two
Intel(R) Xeon(R) Silver 4214 processors were performed. Simulated proton-
proton collisions were used in order to mimic the situation of real data
taking. Depending on the configuration of the outlier removal strategy, an
overall speedup factor between 4 and 5.5 with respect to the default Kalman
filter was achieved. The largest speedup is achieved when no iterations for
the outlier removal are performed. Singling out the calculation steps of the
Kalman filter, i.e. neglecting the part of the algorithms where the
measurement information is constructed, the speedup factor is even larger and
ranges from 5.7 to 10.

In the case of the parametrized Kalman filter, and singling out again the
calculation step of the Kalman filter, $50\%$ of the time is spent
extrapolating the states between the detector layers. Here, the extrapolation
between the UT and the SciFi constitutes the biggest component with a relative
fraction of $40\%$. The remaining Kalman filter steps, consisting of updating
the states with the cluster information and the combination of upstream and
downstream filtered states, are responsible for $16\%$ and $14\%$ of the time
spent, respectively. The extrapolation to the beam line, which is based on the
default LHCb extrapolation algorithm, is responsible for the remaining $20\%$
of the time budget.

##  6 Conclusion

We presented an alternative implementation of a Kalman filter for the LHCb
experiment. Based on simple parametrizations of material effects and the
extrapolation through the magnetic field of the detector, this algorithm
achieves a significant speedup with respect to the current implementation,
while retaining comparable quality of the track parameters. In the future,
further improvements of the parametrizations might allow an even better
estimate of the track parameters and a subsequent speedup. Ideas currently
under discussion include for example an analytic parametrization of the $x$
and $y$ dependence of the parameters employed in the extrapolation from the UT
to the SciFi detector and a better account for the limited acceptance of low
momentum particles. The version presented in this document or a future
implementation might therefore be well suited for the usage in the LHCb
software trigger system for Run 3 of the LHC.

## Acknowledgements

The authors would like to thank the LHCb computing and simulation teams for
their support and for producing the simulated LHCb samples used in the paper.
We also would like to thank the LHCb RTA team for supporting this publication
and reviewing the work. M. De Cian acknowledges support from the Swiss
National Science Foundation grant “Probing right-handed currents in quark
flavour physics”, PZ00P2_174016.

## References

  * [1] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [2] R. Aaij et al., _Performance of the LHCb trigger and full real-time reconstruction in Run 2 of the LHC_ , JINST 14 (2019) P04013, arXiv:1812.10790
  * [3] S. Borghi, _Novel real-time alignment and calibration of the lhcb detector and its performance_ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 845 (2017) 560 , Proceedings of the Vienna Conference on Instrumentation 2016
  * [4] LHCb collaboration, _Framework TDR for the LHCb Upgrade: Technical Design Report_ , CERN-LHCC-2012-007, 2012
  * [5] LHCb collaboration, _LHCb VELO Upgrade Technical Design Report_ , CERN-LHCC-2013-021, 2013
  * [6] LHCb collaboration, _LHCb PID Upgrade Technical Design Report_ , CERN-LHCC-2013-022, 2013
  * [7] LHCb collaboration, _LHCb Tracker Upgrade Technical Design Report_ , CERN-LHCC-2014-001, 2014
  * [8] LHCb collaboration, _LHCb Trigger and Online Technical Design Report_ , CERN-LHCC-2014-016, 2014
  * [9] R. Kalman, _A new approach to linear filtering and prediction problems_ , Journal of Basic Engineering 35 (1960)
  * [10] R. Frühwirth, _Application of kalman filtering to track and vertex fitting_ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 262 (1987) 444
  * [11] E. Bos and E. Rodrigues, _The LHCb Track Extrapolator Tools_ , Tech. Rep. LHCb-2007-140. CERN-LHCb-2007-140, CERN, Geneva, 2007
  * [12] LHCb collaboration, _LHCb magnet: Technical Design Report_ , CERN-LHCC-2000-007, 2000
  * [13] LHCb collaboration, _LHCb RICH: Technical Design Report_ , CERN-LHCC-2000-037, 2000
  * [14] LHCb collaboration, _LHCb calorimeters: Technical Design Report_ , CERN-LHCC-2000-036, 2000
  * [15] LHCb collaboration, _LHCb muon system: Technical Design Report_ , CERN-LHCC-2001-010, 2001
  * [16] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [17] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [18] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [19] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [20] P. Golonka and Z. Was, _PHOTOS Monte Carlo: A precision tool for QED corrections in $Z$ and $W$ decays_, Eur. Phys. J. C45 (2006) 97, arXiv:hep-ph/0506026
  * [21] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [22] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [23] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023