# Particle creation in left-handed metamaterial transmission lines

Alessandro Ferreri  ID Institute for Quantum Computing Analytics (PGI-12),
Forschungszentrum Jülich, 52425 Jülich, Germany  Frank K. Wilhelm  ID
Institute for Quantum Computing Analytics (PGI-12), Forschungszentrum Jülich,
52425 Jülich, Germany

###### Abstract

Transmission lines are excellent examples of quantum simulators of quantum
fields. By properly driving specific circuit elements, these devices can
reproduce quantum relativistic effects such as the particle creation due to
the non-adiabatic stimulation of the quantum vacuum. In this letter, we
investigate the particle creation in left-handed metamaterial transmission
lines (LHTLs). Our results show that, due to the peculiar dispersion
relations, the particle production in LHTLs occurs with much more favorable
conditions with respect to the usual right-handed transmission lines (RHTL).

_Introduction._ — Quantum field theory is the robust mathematical apparatus
that underpins our comprehension of quantum relativistic phenomena [1]. Among
its most successful applications in flat and curved spacetime [2, 3, 4],
quantum field theory predicts highly energetic processes such as the dynamical
Casimir effect (DCE) [5, 6, 7, 8], the Unruh effect [9, 10, 11, 12] and the
Hawking radiation [13, 14, 15]. Unfortunately, due to the considerable energy
scales, the direct observation of these processes in natural contexts has not
been accomplished yet. Nevertheless, in the last decade the impressive
advancements in quantum technologies based on quantum simulation platforms
have been leading to the successful fabrication of devices which well mimic
the main features of such highly energetic phenomena [16, 17, 18, 19, 20].

Quantum simulators are powerful tools for the study of quantum processes whose
reproduction and control in laboratory is, in many case, unfeasible [21, 22].
Among the vast range of possible quantum simulating devices, in this letter we
want to focus our attention on one specific class of them, which can well
describe the dynamics of quantum fields in non-adiabatic scenarios, namely
superconducting circuits based on transmission lines (TLs) [23, 24, 25, 26,
27, 28]. In the frameworks of quantum field theory and cosmology, these
platforms find a large number of theoretical and experimental applications in
the engineering of spacetime analogs [29, 30, 20] and the simulation of
phenomena of particle creation [31, 32, 33], such as the DCE [34, 35, 36] and
the Hawing radiation [37, 38, 39].

In general, the dynamics of the quantum scalar fields and the particle
creation is always described in an environments characterized by positive
dielectric constant and magnetic permeability. Yet, in the last decades a
particular devotion has been dedicated in theoretically and experimentally
engineering metamaterials characterized by left-handed dispersion relations
[40, 41, 42]. In left-handed metamaterials, the group velocity has opposite
sign with respect to the Poynting vector, and this stems from the fact that
both the dielectric constant and the magnetic permeability of the medium are
negative [43].

In one-dimensional systems such as the TL, the left-handedness emerges when
the dispersion law displays a negative group velocity. A standard left-handed
transmission line (LHTL), as that depicted in Fig.1a, can be assembled
starting from the model of a RHTL by interchanging capacitors and inductors
[44] (see also Fig.1a in Sup. Mat.). Due to the mathematical and physical
consequences of this trade, LTHLs can be seen as the dual of the RHTLs [45,
46]. Among other uses, LTHLs and hybrid platforms can find applications in
circuit QED for the simulations of multimode quantum systems [46, 47, 48, 49].

_Theoretical background_ — In this letter, we propose two platforms based on
left-handed transmission lines, whose dynamics is regulated by a set of
superconducting quantum interference devices (SQUIDs) [50, 51]. The difference
between the two schemes relies on the arrangement of the SQUIDs. This implies
two drastically different spectra. We will explore the dynamics of such TLs by
studying the quantum magnetic flux $\Phi_{n}$ (details about the protocol for
the quantization of the TLs are reported in the Sup. Mat.), showing that we
can stimulate the quantum vacuum of the transmission line and generate boson
pairs. Interestingly, we can attribute a different physical origin to the
particle production depending on the arrangement of the circuit elements.

In what follows, we will consider LHTLs made of $N$ nodes. At each node $n$ of
the TL, the SQUID is characterized by the capacitance $C_{\textrm{J}}$, the
Josephson energy $E(t)$, and the phase
$\varphi_{n}=2\pi\Phi_{\textrm{J}n}/\phi_{0}$, where $\Phi_{\textrm{J}n}$ is
the magnetic flux at the SQUID and $\phi_{0}=\pi\hbar/e$, is the magnetic flux
quantum. Importantly, we will assume small amplitude of the plasma oscillation
in the SQUID, $\Phi_{\textrm{J}}/\phi_{0}\ll 1$, and that all SQUIDs work in
the phase regime $E(t)\gg(2e)^{2}/2C_{\textrm{J}}$ [35, 52], thereby expanding
the Josephson energy at the lowest order in $\varphi_{n}/\phi_{0}$ [53].
Finally, we make the identification $\Phi_{\textrm{J}n}=\Phi_{n}$ [36, 52].

The Josephson energy can be externally tuned in order to induce a time-
dependence of dispersion relation. In each platforms discussed in this work,
we will always assume to modulate the Josephson energy around the constant
value $E_{0}=I_{\textrm{c}}\phi_{0}$ via $E(t)=E_{0}[1+2\eta\sin(\Omega
t)]^{2}$ [5], where $I_{\textrm{c}}$ is critical current, $\Omega$ is the
oscillation frequency and $\eta\ll 1$ is dimensionless oscillation amplitude.
To distinguish the three LHTLs we will always refer to Fig.1a as standard
scheme, to Fig.1b as circuit 1, and to Fig.1c as circuit 2.

_Circuit 1: SQUID in parallel._ — As first, we want to study the left-handed
transmission line pictorially represented in Fig.1b. At each node we replace
the inductors of the standard LHTL shown in Fig.1a with a SQUID having
capacitance and Josephson energy $C_{\textrm{J}}$ and $E_{\textrm{J}}=E(t)$
respectively, whereas a capacitor with capacitance $C$ is placed between two
nodes.

Figure 1: Schematic representations of the left-handed transmission line. (a)
Standard LHTL. (b) Circuit 1: LHTL with a set of SQUIDs in parallel. (c)
Circuit 2: LHTL with a set of SQUIDs in series.

In the linearized limit of the SQUID, the Lagrangian of the system expressed
in terms of the magnetic flux is

$\displaystyle\mathcal{L}=\frac{1}{2}\sum_{n=1}^{N}\left[C(\dot{\Phi}_{n}-\dot{\Phi}_{n+1})^{2}+C_{\textrm{J}}\dot{\Phi}_{n}^{2}-\tilde{E}(t)\Phi_{n}^{2}\right],$
(1)

where $\Phi_{n}$ indicates the magnetic flux field at each node, and
$\tilde{E}(t)=(2\pi/\phi_{0})^{2}E(t)$.

In order to describe the dynamics of the transmission line, we need to solve
the equations of motion, namely the Euler-Lagrange equations
$(\textrm{d}/\textrm{d}t)\partial\mathcal{L}/\partial\dot{\Phi}_{n}=\partial\mathcal{L}/\partial\Phi_{n}$
(see the Sup. Mat.). These equations admit a set of solutions of the form
$e^{in\Delta x\,k-i\omega t}$. In order to find the dispersion relation of the
transmission line, we substitute this solution into the equation of motion,
obtaining

$\omega_{j}(t)=\sqrt{\frac{\tilde{E}(t)}{4C\sin^{2}\left(\frac{k_{j}\Delta
x}{2}\right)+C_{\textrm{J}}}},$ (2)

where the wave vector, $k_{j}=\frac{2\pi j}{N\Delta x}$, is bounded within the
first Brillouin zone $j=1,...,N/2$ [46]. The “left-handedness” of the TL
clearly emerges from the group velocity defined as
$v_{\textrm{g}}=\partial\omega/\partial k$, which is negative. Moreover, the
infrared limit of the TL is
$\omega_{\textrm{IR}}=\sqrt{\tilde{E}/(4C+C_{\textrm{J}})}$, which is reached
at the border of the first Brillouin zone, $k_{j}=\pi/\Delta x$.

The modulation of the Josephson energy can simulate the oscillation of the
effective length of the transmission line. To prove it, we assume that the
capacitance of the Josephson junctions is negligible, $C\gg
C_{\textrm{J}}\simeq 0$, moreover, we momentary focus on long wavelength
modes, $k_{j}\ll 1/\Delta x$, and expand the sin function in the denominator
of Eq.(2). This yields

$\displaystyle\omega_{j}(t)\approx$
$\displaystyle\frac{\sqrt{\tilde{E}(t)}}{k_{j}\Delta
x\sqrt{C}}=\frac{N\sqrt{\tilde{E}(t)}}{2\pi
j\sqrt{C}}=\frac{N(t)\sqrt{\tilde{E}_{0}}}{2\pi j\sqrt{C}},$ (3)

where in the last equivalence we moved the modulation from the Josephson
energy to the number of nodes, $N(t)=N\lvert 1+2\eta\sin(\Omega t)\rvert$,
determining an effective length modulation of the LHTL [7, 5].

We note that, depending on the degree of accuracy of our simulation, the
analogy between time-dependent Josephson energy and the modulation of the
effective length could not be accurate. Indeed, although the result in Eq.(3)
recalls the variation in time of the effective length of the TL, the
modulation of $N$ appears only in the dispersion relation, but not in all wave
vectors $k_{j}$ of the mode expansion, as expected by a real modulation of the
cavity length in equivalent cavity QED systems [54, 55]. For this reason, the
TL is in fact simulating a scenario wherein the time-dependence concerns the
speed of light, rather than the cavity length [56, 33]. Nevertheless, although
the interpretation of the modulation of the Josephson energy as time-dependent
speed of light is mathematically more accurate (and can be exploited in
simulations of quantum field theory phenomena in curved spacetime [2]), we
will see that number of bosons generated at the end of the dynamics shows very
strong similarities to the photon number expected in the DCE [5].

_Circuit 2: SQUID in series._ — The second scheme we want to analyze is
illustrated in Fig.1c. In this scheme, we replaced each capacitor of the
standard LHTL in Fig.1a with a SQUID having capacitance $C_{\textrm{J}}$ and
Josephson energy $E(t)$. The Lagrangian of the TL in the linearized limit of
the SQUID is

$\displaystyle\mathcal{L}=$
$\displaystyle\frac{1}{2}\sum_{n=1}^{N}\left[C(\dot{\Phi}_{n+1}-\dot{\Phi}_{n})^{2}-\frac{\Phi_{n}^{2}}{L}-\tilde{E}(t)(\Phi_{n+1}-\Phi_{n})^{2}\right],$
(4)

where we made the identification $C\equiv C_{\textrm{J}}$, and $L$ is the
inductance. As before, assuming a solution of the form $e^{in\Delta
x\,k-i\omega t}$, the Euler-Lagrange equations yield the dispersion relation

$\omega_{j}(t)=\sqrt{\frac{1}{4CL\sin^{2}\left(\frac{k_{j}\Delta
x}{2}\right)}+\frac{\tilde{E}(t)}{C}},$ (5)

with wave vector $k_{j}=\frac{2\pi j}{N\Delta x}$.

Note that this dispersion relation structurally differs from Eq.(2), and in
fact it describes a different left-handed quantum field. Indeed, the
frequencies in Eq.(5) consists of two parts; the first term inside the square
root sets the standard relation between the frequency and wave vector in
LHTLs; while the second term does not depend on the wave vector, and it is
identical for each mode. In the language of quantum field theory, the presence
of latter is the signature of a massive field [1]. In this scheme therefore,
quantum excitations of the magnetic flux field are massive particles with a
left-handed group velocity and quadratic mass proportional to
$M^{2}\propto\tilde{E}(t)/C$. Note that, in contrast to the standard (right-
handed) dispersion relation of the Klein-Gordon field wherein the massive term
of the field becomes negligible in the ultrarelativistic limit ($k\gg
Mc/\hbar$, with $c$ speed of light), in this LHTL the massive term becomes
negligible in long wavelength modes, $k_{j}^{2}\ll 1/[(\Delta
x)^{2}\,CL\tilde{E}(t)]$.

_Frequency and energy of the modes._ — One of the most interesting facts
concerning LHTLs is that the rapid enhancement of the frequency in long
wavelength modes does not reflect the trend of the energy spectrum. Indeed,
regardless of the scheme we are interested in, the Hamiltonian of the LHTL
reads
$\mathcal{\hat{H}}=\sum_{j}\epsilon_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}$,
where the annihilation and creation operators $\hat{a}_{j}$ and
$\hat{a}_{j}^{\dagger}$ fulfill the standard commutation relations
$[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=\delta_{ij}$, and the Hamiltonian
eigenvalues are $\epsilon_{j}=4\hbar\omega_{j}\sin^{2}\left(k_{j}\Delta
x/2\right)$ (see Eqs.(10) and (11) in Sup. Mat.). Due to this discrepancy
between eigenenergies and frequencies in LHTLs, the eigenenergies of LHTLs and
RHTLs behave in a similar manner: the lower the wave vectors, the lower the
eigenenergies. This is evident in Fig.2, where we plot the eigenenergies
$\epsilon_{j}$ of the LHTLs for both circuit 1 (blue), and circuit 2 (cyan),
as well as the dispersion relation in Eq.(34) (red) and Eq.(35) (orange) of
the RHTLs illustrated respectively in Fig.1b and Fig.1c of the Supplemental
Material, together with the dispersion relations in Eq.(2) (purple) and Eq.(5)
(black). We remind that, in contrast with LHTLs, the eigenenergies of the
Hamiltonian in RHTLs coincide with the mode frequencies. This graph shows the
discrepancy between frequencies and eigenenergies in LHTLs: at lower mode
number the eigenenergies tend to smaller values (similarly to the RHTL),
whereas the frequencies increase dramatically.

Figure 2: Hamiltonian eigenenergies of the LHTLs and the RHTLs, and dispersion
relations of the LHTLs, in circuit 1 and 2. Plots refer to the eigenenergies
of the LHTL in circuit 1 (blue), and in circuit 2 (cyan), eigenenergies
(frequencies) of the RHTL in circuit 1 (red), and in circuit 2 (orange),
frequencues of the LHTL in circuit 1 (black), and in circuit 2 (purple).
Indices (1) and (2) label the scheme the curve refers to. Parameters are:
$C=0.4$ pF, $L=60$ pF, $I_{\textrm{c}}=1.25$ $\mu$A and $N=200$.

_Particle creation._ — The time-dependent Josephson energy acts as an external
drive in the dynamics of the transmission lines. When the modulation frequency
$\Omega$ is of the same order of magnitude of the TL mode frequencies,
phenomena of particle creation can occur [2]. To describe such phenomena, we
need to solve the equation of motion with respect to the magnetic flux field
modes, or in other words, we need to find the relation between the
annihilation and creation operators of the magnetic flux modes before the
beginning of the dynamics (input operators) and the ladder operators emerging
once the modulation of the Josephson energy ceases (output operators). Since
the Hamiltonian describes a linear dynamics for small modulations of the
Josephson energy, the output ladder operators are related to the input ones by
Bogoliubov transformations of the form [2, 3, 4]

$\displaystyle\hat{a}_{i}^{\textrm{out}}=\sum_{j}\left(\hat{a}_{j}^{\textrm{in}}\alpha_{ji}+\hat{a}_{j}^{{\dagger}\textrm{in}}\beta_{ji}^{*}\right),$
(6)

where the explicit expression of the coefficients $\alpha_{jh}$ and
$\beta_{jh}$ depend on the specific TL. A more detailed description of the
input-out formalism, as well as the solution of the equations of motion
achieved by means of the multiscale analysis, is reported in the Sup. Mat.

The phenomenon of particle creation is strongly connected to the coefficients
$\beta_{jh}$ in Eq.(16). Indeed, when these coefficients are non-zero, the
initial vacuum state of the quantum field does not correspond to the quantum
vacuum at the end of the dynamics, and the output operators in Eq.(16) do not
act as annihilation operators of the initial quantum vacuum [2]. In the TLs
under considerations, when the modulation frequency of the Josephson energy is
exactly twice the frequency of one of the TL modes, $\Omega=2\omega_{h}$, the
coefficient $\beta_{hh}^{*}$ in Eq.(16) does not vanish, and the operator
$\hat{a}_{h}^{\textrm{out}}$ acts on the input state as a squeezed
annihilation operator. If the system is initially prepared in the vacuum
state, $\lvert 0\rangle=\lvert 0_{1},0_{2},...0_{N/2}\rangle$, we can estimate
the number of particles in the resonant mode $\omega_{h}$ created during the
squeezing process as the expectation value of the output number operator
$\hat{N}_{h}=\hat{a}_{h}^{{\dagger}\textrm{out}}\hat{a}_{h}^{\textrm{out}}$.
Therefore, the output number of particles $\langle\hat{N}_{h}(\tau)\rangle$
and the average energy $\langle\mathcal{\hat{H}}(\tau)\rangle$ at time
$\tau=\eta t$ are respectively given by
$N_{h}(\tau)=\sum_{j}\lvert\beta_{jh}(\tau)\rvert^{2}=\lvert\beta_{hh}(\tau)\rvert^{2}$
and $E_{h}(\tau)=\epsilon_{h}N_{h}(\tau)$.

Before studying these quantities numerically, we first notice that the
analytical formulas describing the output photon number $N_{h}(\tau)$ do not
vary according to the left- or right-handedness of the TL, but it only depends
on the massive or massless behavior of the Josephson energy in the dispersion
relation (see Sup. Mat.). In particular, the number of particles created in a
massless field is given by

$\displaystyle N_{h}(\tau)=\sinh^{2}\left(\kappa_{0h}\tau\right),$ (7)

whereas it corresponds to

$\displaystyle
N_{h}(\tau)=\sinh^{2}\left(\frac{\tilde{E}_{0}\tau}{C\kappa_{0h}}\right)$ (8)

when the Josephson energy defines the mass of the quantum field, with
$\kappa=\omega,\upsilon$. Interestingly, Eq.(7) is formally identical to the
number of photons created via DCE in a cavity confining a quantum scalar field
[5].

The result in Eq.(8) can be interpreted in the key of quantum field theory as
the particle creation due to the modulation of the mass of the quantum field.
In right-handed quantum fields, the time modulation of the massive term can be
associated to the time-dependence of the metric describing the structure of
the spacetime [2, 31, 32]. In case of the LTHL in circuit 2, Eq.(4) does not
correspond to the Lagrangian of a Klein-Gordon scalar field, and the
investigation of gravitational effects on quantum fields with left-handed
dispersion relation would require a more comprehensive analysis which is far
beyond the scope of this letter.

We now want to analyze both $N_{h}(\tau)$ and $E_{h}(\tau)$ in the two LTHLs
depicted in Fig.1b and Fig.1c, as well as in the two RHTLs in Fig.1b and
Fig.1c of the Sup. Mat. Results of our investigation are plotted in Fig.3. In
these graphs, each value represents the expected particle number and the
energy of the mode $h$ at time $\tau$, assuming that this mode is resonant
with the Josephson energy of the SQUID via the resonance condition
$\Omega=2\omega_{h}$ (or $\Omega=2\upsilon_{h}$ for the RHTLs, with
$\upsilon_{h}$ mode frequency).

Figure 3: (a) Expected population and (b) energy at time $\tau$ of the mode
$j$ by activating the resonance condition $\Omega=2\omega_{j}$ in the LHTL and
$\Omega=2\upsilon_{j}$ in the RHTL. Plots refer to LHTL in circuit 1 (blue),
LHTL in circuit 2 (cyan), RHTL in circuit 1 (red), RHTL in circuit 2 (orange).
Indices (1) and (2) label the scheme the curve refers to. Parameters are:
$C=0.4$ pF, $L=60$ pF, $I_{\textrm{c}}=1.25$ $\mu$A, $N=200$, and $\tau=1$ps.

The graph in Fig.3a shows the different trend of the particle creation in
LHTLs and RHTLs. We observe that, in case of massless excitations, the
resonant mode $\upsilon_{h}$ of the RHTL averagely emits more particles at
higher mode numbers (orange dots), while in LHTL the number of massless
particles in the resonant mode $\omega_{h}$ drastically increases at lower
mode numbers (blue dots). This is due to the fact that the function
$N_{h}(\tau)$ in Eq.(7) enhances with the mode frequency (and not with the
mode eigenenergy). The same graph shows that the creation of highly energetic
massive particles is higher in LHTLs (cyan dots) than in RHTLs (red dots).
These results have remarkable consequences in the simulation of particle
creation processes. Indeed, this means that in LHTLs the creation of massless
particles with low energetic modes, and the creation of massive particles with
highly energetic modes are strongly favorable. The particle creation with such
mode numbers by using standard RHTLs is, in contrast, much more difficult to
achieve.

The peculiar trend of the boson number plotted in Fig.3a has also important
consequence in the energy released by the system. Although the eigenvalues of
the Hamiltonian in both LTHLs and RHTLs displayed strong similarities (see
Fig.2), Fig.3b shows that, at lower mode numbers, the left-handed circuit 1
releases a higher output energy with respect to the other platforms. In
particular, this circuit well converts the Josephson energy into quanta of the
magnetic flux, despite the small eigenenergy of these excitations.
Consequently, we can exploit this platform to extract the energy stored in the
Josephson junctions of the SQUIDs with higher efficiency.

_Conclusion._ — This work aimed at studying the quantum dynamics of left-
handed metamaterial transmission lines. Our results demonstrate not only that
the particle production in LTHLs and RHTLs shows strong mathematical
similarities, but also that, due to the peculiar dispersion relation, this
phenomenon is drastically amplified in LHTLs at lower energetic (and therefore
more accessible) modes. These platforms could therefore find interesting
applications in sensing and amplifying low-frequency signals. This work also
paves the way to future investigations of quantum field simulators based on
left-handed metamaterial transmission lines.

_Acknowledgments._ — The authors thank David Edward Bruschi for the helpful
comments. The authors acknowledge support from the joint project No. 13N15685
“German Quantum Computer based on Superconducting Qubits (GeQCoS)” sponsored
by the German Federal Ministry of Education and Research (BMBF) under the
framework “Quantum technologies–from basic research to the market”.

## References

  * Schwartz [2014] M. D. Schwartz,  _Quantum field theory and the standard model_ (Cambridge university press, 2014).
  * Birrell and Davies [1984] N. D. Birrell and P. C. W. Davies,  _Quantum fields in curved space_ (Cambridge university press, 1984).
  * Ford [2021] L. H. Ford, Reports on Progress in Physics 84, 116901 (2021).
  * Martín-Martínez and Menicucci [2014] E. Martín-Martínez and N. C. Menicucci, Classical and Quantum Gravity 31, 214001 (2014).
  * Dodonov [2020] V. Dodonov, Physics 2, 67 (2020).
  * Law [1994] C. K. Law, Physical Review Letters 73, 1931 (1994).
  * Dodonov and Klimov [1996] V. V. Dodonov and A. B. Klimov, Physical Review A 53, 2664 (1996).
  * Crocce  _et al._ [2001] M. Crocce, D. A. R. Dalvit, and F. D. Mazzitelli, Physical Review A 64, 013808 (2001).
  * Unruh [1976] W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * Crispino  _et al._ [2008] L. C. B. Crispino, A. Higuchi, and G. E. A. Matsas, Rev. Mod. Phys. 80, 787 (2008).
  * Fewster  _et al._ [2016] C. J. Fewster, B. A. Juárez-Aubry, and J. Louko, Classical and Quantum Gravity 33, 165003 (2016).
  * Bruschi  _et al._ [2010] D. E. Bruschi, J. Louko, E. Martín-Martínez, A. Dragan, and I. Fuentes, Phys. Rev. A 82, 042332 (2010).
  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Parikh and Wilczek [2000] M. K. Parikh and F. Wilczek, Phys. Rev. Lett. 85, 5042 (2000).
  * Almheiri  _et al._ [2021] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini, Rev. Mod. Phys. 93, 035002 (2021).
  * Wilson  _et al._ [2010] C. M. Wilson, T. Duty, M. Sandberg, F. Persson, V. Shumeiko, and P. Delsing, Physical Review Letters 105, 233907 (2010).
  * Wilson  _et al._ [2011] C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, and P. Delsing, Nature 479, 376 (2011).
  * Steinhauer [2016] J. Steinhauer, Nature Physics 12, 959 (2016).
  * Muñoz De Nova  _et al._ [2019] J. R. Muñoz De Nova, K. Golubkov, V. I. Kolobov, and J. Steinhauer, Nature 569, 688 (2019).
  * Shi  _et al._ [2023] Y.-H. Shi, R.-Q. Yang, Z. Xiang, Z.-Y. Ge, H. Li, Y.-Y. Wang, K. Huang, Y. Tian, X. Song, D. Zheng, K. Xu, R.-G. Cai, and H. Fan, Nature Communications 14, 3263 (2023).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * You and Nori [2011] J. Q. You and F. Nori, Nature 474, 589 (2011).
  * Yurke and Denker [1984] B. Yurke and J. S. Denker, Physical Review A 29, 1419 (1984).
  * Blais  _et al._ [2004] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Physical Review A 69, 062320 (2004).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. Girvin, and A. Wallraff, Reviews of Modern Physics 93, 025005 (2021).
  * Zhang  _et al._ [2023] J. Zhang, R. Ferguson, S. Kühn, J. F. Haase, C. Wilson, K. Jansen, and C. A. Muschik, Quantum 7, 1148 (2023).
  * Pozar [2011] D. M. Pozar,  _Microwave engineering_ (John wiley & sons, 2011).
  * Vool and Devoret [2017] U. Vool and M. Devoret, International Journal of Circuit Theory and Applications 45, 897 (2017).
  * Sabín [2022] C. Sabín, Universe 8, 10.3390/universe8090452 (2022).
  * Terrones and Sabín [2021] A. Terrones and C. Sabín, Universe 7, 499 (2021).
  * Nation  _et al._ [2012] P. D. Nation, J. R. Johansson, M. P. Blencowe, and F. Nori, Reviews of Modern Physics 84, 1 (2012).
  * Tian  _et al._ [2017] Z. Tian, J. Jing, and A. Dragan, Physical Review D 95, 125003 (2017).
  * Lang and Schützhold [2019] S. Lang and R. Schützhold, Physical Review D 100, 065003 (2019).
  * Lähteenmäki  _et al._ [2013] P. Lähteenmäki, G. S. Paraoanu, J. Hassel, and P. J. Hakonen, Proceedings of the National Academy of Sciences 110, 4234 (2013).
  * Johansson  _et al._ [2009] J. R. Johansson, G. Johansson, C. M. Wilson, and F. Nori, Physical Review Letters 103, 147003 (2009).
  * Johansson  _et al._ [2010] J. R. Johansson, G. Johansson, C. M. Wilson, and F. Nori, Physical Review A 82, 052509 (2010).
  * Nation  _et al._ [2009] P. D. Nation, M. P. Blencowe, A. J. Rimberg, and E. Buks, Physical Review Letters 103, 087004 (2009).
  * Tian and Du [2019] Z. Tian and J. Du, The European Physical Journal C 79, 994 (2019).
  * Blencowe and Wang [2020] M. P. Blencowe and H. Wang, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378, 20190224 (2020).
  * Smith  _et al._ [2000] D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, Phys. Rev. Lett. 84, 4184 (2000).
  * Smith and Kroll [2000] D. R. Smith and N. Kroll, Phys. Rev. Lett. 85, 2933 (2000).
  * Houck  _et al._ [2003] A. A. Houck, J. B. Brock, and I. L. Chuang, Phys. Rev. Lett. 90, 137401 (2003).
  * Veselago [1967] V. G. Veselago, Usp. fiz. nauk 92, 517 (1967).
  * Kozyrev and Van Der Weide [2008] A. B. Kozyrev and D. W. Van Der Weide, Journal of Physics D: Applied Physics 41, 173001 (2008).
  * Alizadeh  _et al._ [2020] A. Alizadeh, B. Rejaei, and M. Fardmanesh, IEEE Transactions on Applied Superconductivity 30, 1 (2020).
  * Egger and Wilhelm [2013] D. J. Egger and F. K. Wilhelm, Physical Review Letters 111, 163601 (2013).
  * Messinger  _et al._ [2019] A. Messinger, B. G. Taketani, and F. K. Wilhelm, Physical Review A 99, 032325 (2019).
  * Wang  _et al._ [2019] H. Wang, A. Zhuravel, S. Indrajeet, B. Taketani, M. Hutchings, Y. Hao, F. Rouxinol, F. Wilhelm, M. LaHaye, A. Ustinov, and B. Plourde, Physical Review Applied 11, 054062 (2019).
  * Indrajeet  _et al._ [2020] S. Indrajeet, H. Wang, M. Hutchings, B. Taketani, F. K. Wilhelm, M. LaHaye, and B. Plourde, Physical Review Applied 14, 064033 (2020).
  * Ovchinnikova  _et al._ [2013] E. A. Ovchinnikova, S. Butz, P. Jung, V. P. Koshelets, L. V. Filippenko, A. S. Averkin, S. V. Shitov, and A. V. Ustinov, Superconductor Science and Technology 26, 114003 (2013).
  * Shramkova  _et al._ [2017] O. V. Shramkova, N. Lazarides, G. P. Tsironis, and A. V. Ustinov, Applied Physics A 123, 58 (2017).
  * Weißl  _et al._ [2015] T. Weißl, B. Küng, E. Dumur, A. K. Feofanov, I. Matei, C. Naud, O. Buisson, F. W. J. Hekking, and W. Guichard, Physical Review B 92, 104508 (2015).
  * Wallquist  _et al._ [2006] M. Wallquist, V. S. Shumeiko, and G. Wendin, Phys. Rev. B 74, 224506 (2006).
  * Law [1995] C. K. Law, Physical Review A 51, 2537 (1995).
  * Ferreri  _et al._ [2022] A. Ferreri, H. Pfeifer, F. K. Wilhelm, S. Hofferberth, and D. E. Bruschi, Phys. Rev. A 106, 033502 (2022).
  * Sabín [2018] C. Sabín, New Journal of Physics 20, 053028 (2018).
  * Bender  _et al._ [1999] C. M. Bender, S. Orszag, and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory_ , Vol. 1 (Springer Science & Business Media, 1999).

##  I Supplemental Material to Particle creation in left-handed metamaterial
transmission line

###### Contents

  1. I Supplemental Material to Particle creation in left-handed metamaterial transmission line
  2. II Equation of motion, Hamiltonian and quantization
  3. III Dynamics and Bogoliubov transformation
  4. IV Multiple scale analysis
  5. V Right-handed transmission lines

##  II Equation of motion, Hamiltonian and quantization

In this section we want to discuss the formalism employed for the quantization
of the two LHTLs presented in this work, starting from first principles. Our
goal is therefore to achieve a quantized expression of both the modes of the
quantum magnetic flux fields and the Hamiltonian at $t<0$, namely before the
beginning of the modulation of the Josephson energy. For this reason, we
conveniently omit the time-dependence in the Josephson energy, and write
$\tilde{E}=E_{0}(2\pi/\phi_{0})^{2}$.

As first step, we need to solve the equations of motion of the two TLs. These
can be achieved from the two Lagrangians in Eq.(1) and Eq.(4) of the main
text, respectively for the circuit 1 and the circuit 2. At any node $1<n<N$,
the equations of motion become

$\displaystyle
2\ddot{\Phi}_{n}-\ddot{\Phi}_{n+1}-\ddot{\Phi}_{n-1}+\frac{\tilde{E}}{C}\Phi_{n}$
$\displaystyle=0,$ (1)

for circuit 1, and

$\displaystyle
2\ddot{\Phi}_{n}-\ddot{\Phi}_{n+1}-\ddot{\Phi}_{n-1}+\frac{\Phi_{n}}{CL}+\frac{\tilde{E}}{C}(2\Phi_{n}-\Phi_{n+1}-\Phi_{n-1})=0,$
(2)

for circuit 2. Note that in the left part of Eq.(1) we omitted the term
$\frac{C_{\textrm{J}}}{C}\ddot{\Phi}_{n}$, since we assumed $C\gg
C_{\textrm{J}}\simeq 0$. Moreover, in Eq.(2) we made the identification
$C_{\textrm{J}}\equiv C$.

In our description, both transmission lines consist of a finite number of
nodes $N$. When $N$ is large enough, we can ignore the effect at the edges of
the TLs and restrict the solution to periodic boundary conditions [2].
Assuming solution of the form $e^{in\Delta x\,k-i\omega t}$, the quantum
magnetic flux is described by the mode decomposition

$\displaystyle\Phi(n,t)=\sum_{j=1}^{N/2}\phi_{j}(n,t)\hat{a}_{j}+\textrm{h.c.},$
(3)

where the modes are

$\displaystyle\phi_{j}(n,t\leq
0)=\sqrt{\frac{\hbar}{2CN\omega_{0j}}}e^{i(k_{j}n\Delta x-\omega_{0j}t)},$ (4)

and $\omega_{0j}\equiv\omega_{j}(t\leq 0)$ are the frequencies. Note that,
since the minimum distinguishable wavelength is
$\lambda_{\textrm{min}}=2\Delta x$ [46], the sum over all modes runs up to
$N/2$. In Eq.(3) we quantized the amplitude of the field by means of standard
bosonic commutation relation $[\hat{a}_{i},\hat{a}_{j}]=\delta_{ij}$. Note
that the modes are normalized such that

$\displaystyle-\frac{iC}{\hbar}\sum_{n=1}^{N}\left[\phi_{i}(n,t)\frac{\partial\phi_{j}^{*}(n,t)}{\partial
t}-\frac{\partial\phi_{i}^{*}(n,t)}{\partial
t}\phi_{j}(n,t)\right]=\delta_{ij},$ (5)

where we made use of the representation of the Kronecker delta

$\displaystyle\delta_{lh}=\frac{1}{N}\sum_{n=1}^{N}e^{2\pi in(l-h)/N}.$ (6)

The Hamiltonian of the transmission lines is achieved from the Lagrangian by
means of the Legendre transformation,
$\mathcal{H}=\sum_{n}\frac{\partial\mathcal{L}}{\partial\dot{\Phi}_{n}}\dot{\Phi}_{n}-\mathcal{L},$
thereby obtaining

$\displaystyle\mathcal{H}=$
$\displaystyle\frac{1}{2}\sum_{n=1}^{N}\left[\sum_{i=1}^{n}\frac{(P_{i})^{2}}{C}+\frac{(\Phi_{n})^{2}}{L}\right],$
(7)

with conjugated momenta

$\displaystyle P_{n}=$
$\displaystyle\frac{\partial\mathcal{L}}{\partial\dot{\Phi}_{n}}=C\left(2\dot{\Phi}_{n}-\dot{\Phi}_{{n}+1}-\dot{\Phi}_{{n}-1}\right).$
(8)

By substituting the mode decomposition in Eq.(3) and exploiting the
normalization condition in Eq.(5), the Hamiltonian reduces to

$\displaystyle\mathcal{\hat{H}}=$ $\displaystyle
4\hbar\sum_{j=1}^{N/2}\omega_{0j}\sin^{2}\left(\frac{k_{j}\Delta
x}{2}\right)\hat{a}_{j}^{\dagger}\hat{a}_{j}.$ (9)

Note that the form of both the modes in Eq.(4) and the Hamiltonian in Eq.(9)
is valid for both LHTLs under consideration. The difference between the two
schemes relies only on the dispersion relations and the Hamiltonian
eigenenergies. In particular, the explicit form of the eigenenergies in
circuit 1 is

$\epsilon^{(1)}_{j}=2\hbar\sin\left(\frac{k_{j}\Delta
x}{2}\right)\sqrt{\frac{\tilde{E}(t)}{C}},$ (10)

whereas the eigenenergies in circuit 2 are

$\epsilon^{(2)}_{j}=2\hbar\sin\left(\frac{k_{j}\Delta
x}{2}\right)\sqrt{\frac{1}{CL}+\frac{4\tilde{E}(t)\sin^{2}\left(\frac{k_{j}\Delta
x}{2}\right)}{C}}.$ (11)

##  III Dynamics and Bogoliubov transformation

Arguments in this section are valid for both LHTLs considered in this work. To
study the dynamics at $t>0$ we will use the procedure employed in [8, 57].
When the modulation of the Josephson energy starts, each mode can be written
as:

$\displaystyle\phi_{j}(n,t>0)=\sum_{l=1}^{N/2}\sqrt{\frac{\hbar}{CN}}Q_{jl}(t)e^{ik_{l}n\Delta
x},$ (12)

where $Q_{jl}(t)$ are solution of the equations of motion with continuity
conditions $Q_{jl}(0)=\delta_{jl}/\sqrt{2\omega_{0j}}$ and
$\dot{Q}_{jl}(0)=-i\delta_{jl}\sqrt{2\omega_{0j}}$.

The explicit expression for $Q_{jl}(t)$ are achieved by solving the equations
of motion. The substitution of Eq.(12) into the equations of motion yields:

$\displaystyle\sum_{l=1}^{N/2}\left[\frac{\ddot{Q}_{jl}(t)}{\omega_{j}^{2}(t)}+Q_{jl}(t)\right]e^{in\Delta
xk_{l}}$ $\displaystyle=0.$ (13)

By multiply both sides of the equation by $e^{-in\Delta xk_{h}}/N$, summing
over $n$, and exploiting the representation of the Kronecker delta in Eq.(6),
we get

$\displaystyle\ddot{Q}_{jh}(t)+\omega_{h}^{2}(t)Q_{jh}(t)=0.$ (14)

To solve it, we employ the multiple scale analysis. This strategy, as well as
the solutions of this equation, are reported in the next section. Once the
Josephson energy arrests at its original value at time $t_{\textrm{f}}$,
$E(t_{\textrm{f}})=E(0)\equiv E_{0}$, the solution of the equation of motion
simply becomes:

$\displaystyle Q_{jh}=A_{jh}e^{i\omega_{0h}t}+B_{jh}e^{-i\omega_{0h}t},$ (15)

where the coefficients $A_{jh}$ and $B_{jh}$ are determined by the continuity
condition of each $Q_{jh}$ at the end of the motion. The explicit form of the
coefficients $A_{jh}$ and $B_{jh}$ is calculated for the two LHTLs in the next
section.

Note that the $in$ and $output$ modes determine the quantum vacuum at
different time, and their relation is described by Bogoliubov transformations:

$\displaystyle\hat{a}_{h}^{\textrm{out}}=\sum_{j}\left(\hat{a}_{j}^{\textrm{in}}\alpha_{jh}+\hat{a}_{j}^{{\dagger}\textrm{in}}\beta_{jh}^{*}\right).$
(16)

Substituting Eq.(15) into Eq.(12), we achieve the form of output Fourier mode
expansion of the field with respect to the input modes. Recombining all terms
properly and comparing with respect to Eq.(16) we observe that

$\displaystyle\alpha_{jh}=$ $\displaystyle\sqrt{2\omega_{0h}}B_{jh}$
$\displaystyle\beta_{jh}=$ $\displaystyle\sqrt{2\omega_{0h}}A_{jh}.$ (17)

These formulas relate both Bogoliubov coefficients to the coefficients
$A_{jh}$ and $B_{jh}$.

##  IV Multiple scale analysis

In this section, we present the solution strategy to solve Eq.(14) for the two
LHTLs. This formalism is presented in [57]. As first, we define a new “time
scale” $\tau=\eta t$ and expand $Q_{jh}(t)$ with respect to $\eta$:

$\displaystyle Q_{jh}(t)=$ $\displaystyle Q_{jh}^{(0)}(t,\tau)+\eta
Q_{jh}^{(1)}(t,\tau)$ $\displaystyle\ddot{Q}_{jh}(t)=$
$\displaystyle\partial_{t}^{2}Q_{jh}^{(0)}(t,\tau)+\eta\left[2\partial_{\tau
t}^{2}Q_{jh}^{(0)}(t,\tau)+\partial_{t}^{2}Q_{jh}^{(1)}(t,\tau)\right]$ (18)

We now substitute Eq.(IV) into Eq. (14) and solve this equation at different
orders in $\tau$. At the zeroth order in $\tau$ we simply achieve

$\displaystyle\ddot{Q}_{jh}^{(0)}+\omega_{h0}^{2}Q_{jh}^{(0)}=0,$ (19)

whose solution is

$\displaystyle
Q_{jh}^{(0)}(t,\tau)=A_{jh}(\tau)e^{i\omega_{h0}t}+B_{jh}(\tau)e^{-i\omega_{h0}t},$
(20)

and from the continuity conditions for $Q_{jh}(t)$ we have the following
initial conditions for $A_{jh}(\tau)$ and $B_{jh}(\tau)$:

$\displaystyle A_{jh}(0)=$ $\displaystyle 0,$ $\displaystyle B_{jh}(0)=$
$\displaystyle\frac{1}{\sqrt{2\omega_{h0}}}\delta_{jh}.$ (21)

We now focus on the circuit 1, describing a quantum magnetic flux field with
massless excitations. We remind that the dispersion relation for such field is
given by Eq.(2) of the main text at $t=0$, with $C_{\textrm{J}}=0$.
Considering only the first order in $\tau$, Eq.(14) becomes:

$\displaystyle 2\partial_{\tau
t}Q_{jh}^{(0)}+\partial_{t}^{2}Q_{jh}^{(1)}+\omega_{h0}^{2}Q_{jh}^{(1)}+4\omega_{h0}^{2}\sin(\Omega
t)Q_{jh}^{(0)}=0,$ (22)

Substituting the solution of the zeroth order, Eq.(20), into Eq.(22) we obtain

$\displaystyle\partial_{t}^{2}Q_{jh}^{(1)}+\omega_{h0}^{2}Q_{jh}^{(1)}=$
$\displaystyle-2\partial_{\tau t}Q_{jh}^{(0)}-4\omega_{h0}^{2}\sin(\Omega
t)Q_{jh}^{(0)}$ $\displaystyle=$
$\displaystyle-2i\omega_{h0}\left(\frac{\partial
A_{jh}}{\partial\tau}e^{i\omega_{h0}t}-\frac{\partial
B_{jh}}{\partial\tau}e^{-i\omega_{h0}t}\right)+2i\omega_{h0}^{2}\left(e^{i\Omega
t}-e^{-i\Omega
t}\right)\left(A_{jh}e^{i\omega_{h0}t}+B_{jh}e^{-i\omega_{h0}t}\right)$
$\displaystyle=$
$\displaystyle-2i\omega_{h0}e^{i\omega_{h0}t}\left(\frac{\partial
A_{jh}}{\partial\tau}-\omega_{h0}B_{jh}e^{i(\Omega-2\omega_{h0})t}\right)+2i\omega_{h0}e^{-i\omega_{h0}t}\left(\frac{\partial
B_{jh}}{\partial\tau}-\omega_{h0}A_{jh}e^{-i(\Omega-2\omega_{h0})t}\right).$
(23)

We now want the solution of such equation to be lacking of secularities.
Secularities are all those terms proportional to $e^{\pm i\omega_{h0}t}$,
namely all those term which are already solution of the homogeneous equation.
To avoid such singularities, we need the coefficients of $e^{\pm
i\omega_{h0}t}$ to vanish, thereby obtaining the following set of differential
equations:

$\displaystyle\frac{\partial
A_{jh}}{\partial\tau}-\omega_{h0}\delta(\Omega-2\omega_{h0})B_{jh}=$
$\displaystyle 0,$ (24) $\displaystyle\frac{\partial
B_{jh}}{\partial\tau}-\omega_{h0}\delta(\Omega-2\omega_{h0})A_{jh}=$
$\displaystyle 0,$ (25)

or

$\displaystyle\frac{\partial A_{jh}}{\partial\tau}=$
$\displaystyle\omega_{h0}\delta(\Omega-2\omega_{h0})B_{jh},$ (26)
$\displaystyle\frac{\partial B_{jh}}{\partial\tau}=$
$\displaystyle\omega_{h0}\delta(\Omega-2\omega_{h0})A_{jh},$ (27)

Such equations are solved by differentiating one, substituting it into the
other one and exploiting the initial conditions in Eq.(IV). Finally we get:

$\displaystyle A_{jh}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{h0}}}\sinh\left(\omega_{h0}\tau\right)\delta_{jh}\delta(\Omega-2\omega_{h0}),$
$\displaystyle B_{jh}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{h0}}}\cosh\left(\omega_{h0}\tau\right)\delta_{jh}\delta(\Omega-2\omega_{h0}).$
(28)

We proceed similarly for circuit 2, describing a quantum field with massive
excitations. At the zeroth order in $\tau$, the equation for $Q_{jh}$ is
identical to Eq.(19), where now $\omega_{h0}$ corresponds to Eq.(5) of the
main text at $t=0$. At the first order, the equation is:

$\displaystyle 2\partial_{\tau
t}Q_{jh}^{(0)}+\partial_{t}^{2}Q_{jh}^{(1)}+\omega_{h0}^{2}Q_{jh}^{(1)}+\frac{4E_{0}}{C}\sin(\Omega
t)Q_{jh}^{(0)}=0.$ (29)

Moving the zeroth order to the right side and replacing the solution of the
zeroth order, Eq.(20), into Eq.(29) we get

$\displaystyle\partial_{t}^{2}Q_{jh}^{(1)}+\omega_{h0}^{2}Q_{jh}^{(1)}=$
$\displaystyle-2\partial_{\tau
t}Q_{jh}^{(0)}-\frac{4\tilde{E}_{0}}{C}\sin(\Omega t)Q_{jh}^{(0)}$
$\displaystyle=$ $\displaystyle-2i\omega_{h0}\left(\frac{\partial
A_{jh}}{\partial\tau}e^{i\omega_{h0}t}-\frac{\partial
B_{jh}}{\partial\tau}e^{-i\omega_{h0}t}\right)+\frac{2i\tilde{E}_{0}}{C}\left(e^{i\Omega
t}-e^{-i\Omega
t}\right)\left(A_{jh}e^{i\omega_{h0}t}+B_{jh}e^{-i\omega_{h0}t}\right)$
$\displaystyle=$
$\displaystyle-2ie^{i\omega_{h0}t}\left(\omega_{h0}\frac{\partial
A_{jh}}{\partial\tau}-\frac{\tilde{E}_{0}}{C}B_{jh}e^{i(\Omega-2\omega_{h0})t}\right)+2ie^{-i\omega_{h0}t}\left(\omega_{h0}\frac{\partial
B_{jh}}{\partial\tau}-\frac{\tilde{E}_{0}}{C}A_{jh}e^{-i(\Omega-2\omega_{h0})t}\right).$
(30)

As before, to get rid of all secularities we need to solve the couple of
differential equations

$\displaystyle\frac{\partial A_{jh}}{\partial\tau}=$
$\displaystyle\frac{\tilde{E}_{0}}{C\omega_{h0}}\delta(\Omega-2\omega_{h0})B_{jh},$
(31) $\displaystyle\frac{\partial B_{jh}}{\partial\tau}=$
$\displaystyle\frac{\tilde{E}_{0}}{C\omega_{h0}}\delta(\Omega-2\omega_{h0})A_{jh},$
(32)

whose solutions are

$\displaystyle A_{jh}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{h0}}}\sinh\left(\frac{\tilde{E}_{0}\tau}{C\omega_{h0}}\right)\delta_{jh}\delta(\Omega-2\omega_{h0}),$
$\displaystyle B_{jh}=$
$\displaystyle\frac{1}{\sqrt{2\omega_{h0}}}\cosh\left(\frac{\tilde{E}_{0}\tau}{C\omega_{h0}}\right)\delta_{jh}\delta(\Omega-2\omega_{h0}).$
(33)

##  V Right-handed transmission lines

Figure 1: Schematic representations of the right-handed transmission line. (a)
Standard RHTL. (b) Circuit 1: RHTL with a set of SQUIDs in parallel. (c)
Circuit 2: RHTL with a set of SQUIDs in series.

In this section, we want to make a comparison between our two LHTLs and the
equivalent RHTLs illustrated in Fig.1b and Fig.1c. The first scheme we want to
consider, which we will name right-handed circuit 1, is represented in Fig.1b,
and it is achieved from the standard RHTL in Fig.1a by replacing all
capacitors with SQUIDs. A similar circuit was already proposed as analogue of
a scalar field interacting with the the gravitational field [32]. The
dispersion relation of such RHTL is

$\displaystyle\upsilon_{j}(t)=\sqrt{\frac{4\sin^{2}\left(\frac{p_{j}\Delta
x}{2}\right)}{LC}+\frac{\tilde{E}_{\textrm{J}}}{C}},$ (34)

with $C\equiv C_{\textrm{J}}$. In this work, $\upsilon_{j}$ and $p_{j}$ always
indicate the frequency and the wave vector in the RHTL, respectively. In the
framework of quantum field theory, this expression suggests that the second
term inside the square root is related to the presence of a massive term
proportional to $\sqrt{E(t)/C}$ [38].

The second RHTL, which we will refer to as right-handed circuit 2, is
illustrated in Fig.1c and is implemented from the scheme in Fig.1a by
replacing all inductors with a set of SQUIDs. A similar scheme was realized
for the study of nonlinear effects in the bare frequencies of the Hamiltonian
(Kerr effect) [52]. In the regime $C\gg C_{\textrm{J}}\simeq 0$, the
dispersion relation in this TL reads

$\displaystyle\upsilon_{j}(t)=$ $\displaystyle
2\sqrt{\frac{\tilde{E}(t)}{C}}\sin\left(\frac{p_{j}\Delta x}{2}\right).$ (35)

The Hamiltonian of RHTLs under consideration takes the same form for both
schemes in Fig.1b and Fig.1c:

$\displaystyle\mathcal{\hat{H}}_{\textrm{R}}=$
$\displaystyle\hbar\sum_{j=1}^{N/2}\upsilon_{j}\hat{b}_{j}^{\dagger}\hat{b}_{j}.$
(36)

A direct comparison between the dispersion relations in LHTLs (see Eq.(2) and
(5) in the main text), and in RHTLs (see (34) and (35)) shows interesting
analogies. For instance, the Josephson energy can cover the role of a massless
or a massive term in both types of TLs. Finally, the fact that the dispersion
relation of the left-handed circuit 1 shares the same form of the dispersion
relation of the right-handed circuit 2, and analogously for the other two
schemes, is a further signature of the duality between left-handed and right-
handed transmission lines. A crucial consequence of this duality is that, by
repeating the same procedure and solving the equation of motions for RHTLs, we
exactly obtain Eqs.(28) and Eqs.(33) in the massless and massive cases
respectively (with $\upsilon_{0h}$ instead of $\omega_{0h}$).


11institutetext: University of Cincinnati 22institutetext: Princeton
University 33institutetext: Massachusetts Institute of Technology

# Progress in developing a hybrid deep learning algorithm for identifying and
locating primary vertices

Simon Akar 11<EMAIL_ADDRESS>Gowtham Atluri 11 Thomas Boettcher 11 Michael
Peters 11 Henry Schreiner 22 Michael Sokoloff 11<EMAIL_ADDRESS>Marian
Stahl 11 William Tepe 11 Constantin Weisser 33 Mike Williams 33

###### Abstract

The locations of proton-proton collision points in LHC experiments are called
primary vertices (PVs). Preliminary results of a hybrid deep learning
algorithm for identifying and locating these, targeting the Run 3 incarnation
of LHCb, have been described at conferences in 2019 and 2020. In the past year
we have made significant progress in a variety of related areas. Using two
newer Kernel Density Estimators (KDEs) as input feature sets improves the
fidelity of the models, as does using full LHCb simulation rather than the
“toy Monte Carlo” originally (and still) used to develop models. We have also
built a deep learning model to calculate the KDEs from track information.
Connecting a tracks-to-KDE model to a KDE-to-hists model used to find PVs
provides a proof-of-concept that a single deep learning model can use track
information to find PVs with high efficiency and high fidelity. We have
studied a variety of models systematically to understand how variations in
their architectures affect performance. While the studies reported here are
specific to the LHCb geometry and operating conditions, the results suggest
that the same approach could be used by the ATLAS and CMS experiments.

##  1 Introduction

The LHCb experiment is currently being upgraded for the planned start of Run 3
of the LHC in 2022. It will record proton-proton collision data at five times
the instantaneous luminosity of Run 2. The average number of visible primary
vertices (PVs), proton-proton collision points with a least 5 charged tracks
reconstructed in the detector closest to the beam-crossing region, will
increase from 1.1 to 5.6. Building on the success of Run 2
LHCbCollaboration:2014vzo , the experiment will move to a pure software data
ingestion and trigger system, eliminating the Level 0 hardware trigger
altogether. A conventional PV finding algorithm Reiss:2749592 ; Reiss:CtD2020
that satisfies all requirements defined in the Technical Design Report
Aaij:2018jht  serves as the baseline. In parallel, we have been developing a
hybrid machine learning algorithm that may be more performant – providing
better fidelity and executing more quickly.

A cartoon illustrating the Run 3 Vertex Locator (VELO) is shown in Fig. 1.
Approximately 41 million pixels, $55\times 55\,\mu{\rm m}^{2}$ each, will
populate 26 circular discs oriented perpendicular to the beamline. The LHCb
detector is a forward spectrometer with the remaining elements to the right of
the VELO in this view. In the longitudinal direction, the luminous region can
be characterized by the ellipse labeled as the interaction region. In the
transverse directions, most PVs are produced within 40 $\mu$m of the center of
the beamline. An error ellipse with its minor axis drawn to the vertical scale
of the figure would be a factor of 100 thinner than that shown. The typical
longitudinal resolution of a PV reconstructed following full tracking and
fitting is $40-200\,\mu$m (depending primarily on track multiplicity), and the
typical transverse resolution is $\sim 10\,\mu$m. These small dimensions make
it useless to try to visualize them accurately on the scale of the VELO. They
also allow tracking and PV-finding algorithms to use these constraints to
execute quickly and efficiently. Of interest for the work presented here, PV-
finding algorithms can use track parameters evaluated at their points of
closest approach (poca) to the beamline to convert sparse point clouds of
three dimensional data to rich one dimensional data sets amenable to
processing by deep neural networks.

Figure 1: This diagram illustrates the luminous region of the LHCb experiment
and the Vertex Locator (VELO), the high precision silicon pixel detector that
surrounds it. Primary vertices are too small to be seen if drawn to scale. For
details, see the text.

The initial algorithm defined a rich, one-dimensional Kernel Density Estimator
(KDE) histogram, plus two more one-dimensional histograms, to describe the
probabilities of tracks traversing small voxels in space Fang:2019wsd ;
akar2020updated . A convolutional Neural Network (CNN) produces a single one-
dimensional histogram that nominally predicts Gaussian peaks at the locations
of the true PVs using the the three input histograms as its feature set. A
hand-written clustering algorithm then identifies the candidate PVs and their
positions. The results reported earlier used a “toy Monte Carlo” wit proto-
tracking Fang:2019wsd . As discussed below, using track parameters produced by
a proposed LHCb Run 3 Vertex Locator (VELO) tracking algorithm
Hennequin:2019itm  leads to significantly better performance.

The original KDE akar2020updated  is a projection of a three-dimensional
probability distribution in voxels that has contributions only when two tracks
pass close to each other. Calculating this KDE exactly is very time-consuming,
and we want to replace it with a deep learning (DL) algorithm of its own.
Initial attempts suggested that learning track interactions would be
difficult, so we defined two new KDEs, one summing probabilities of individual
tracks in spatial voxels and one summing squares of probabilities. The
original KDE used a combination of these. As discussed below, using only one
of these KDEs as an input feature leads to worse results than using the
original KDE, but using both leads to significantly better performance.

Defining KDEs that depend only on individual track parameters, and not on the
direct interaction of two tracks, allows a tracks-to-KDE DL algorithm to build
a pretty good estimate of the numerically calculated “probability” KDE using
the ensemble of individual tracks parameters as input features. We have
connected a preliminary model of this DL algorithm to a simple, somewhat
limited performance KDE-to-hists model to test whether such a tracks-to-hists
algorithm can achieve as good fidelity as the KDE-to-hists algorithm does
using a numerically calculated KDEs. Initial results are encouraging.

We have studied modified versions of our original CNN architecture
Fang:2019wsd  to understand how performance depends on the number of model
parameters, the number of network layers, the number of channels per layer,
whether batch normalization is used, and the extent to which skip connections
are used. We refer to these models as members of the AllCNN family. We have
also tested a one-dimensional version of the U-Net architecture
ronneberger2015unet , originally developed for two-dimensional biomedical
images, using toy Monte Carlo and find that it (slightly) outperforms the best
of our AllCNN models when both are trained with the same KDE.

##  2 Performance Evolution

Figure 2:  This plot compares the performance of models reported in previous
years (labeled ACAT-2019 and CtD-2020) with new models described in detail in
the text. An asymmetry parameter between the cost of overestimating
contributions to the target histograms and underestimating them is varied to
produce the families of points observed.

Figure 2 shows how the performance of the KDE-to-hists algorithms have evolved
over time. The solid blue circles show the performance of any early model
described at ACAT-2019 Fang:2019wsd . The green squares and magenta diamonds
show the performances described at Connecting-the-Dots in 2020 akar2020updated
. The efficiency is shown on the horizontal axis and the false positive rate
per event is shown on the vertical axis. An asymmetry parameter between the
cost of overestimating contributions to the target histograms and
underestimating them Fang:2019wsd  is varied to produce the families of points
observed. Those models were trained and tested using the original KDE and toy
MC. The reported results come from statistically independent validation
samples. In the fall of 2020 we tested the model of the magenta diamonds using
full LHCb MC in which the (original) KDEs were derived from a full VELO
tracking algorithm Hennequin:2019itm  and we observed that the performance was
better than with toy MC (not shown). We then trained using a full LHCb MC
sample and observed the results plotted as red triangles. We ascribe the clear
improvement in performance (better efficiency for the same false positive
rate, or lower false positive rate for the same efficiency) to a lower “ghost
rate” in the full VELO tracking relative to the proto-tracking in the toy MC
and to better track parameter determinations from using more hits in a more
sophisticated fit. An algorithm that replaces the original KDE with two KDEs
produces even better results, plotted as cyan-filled circles. The difference
is most pronounced at efficiencies greater than $\sim 94\%$. The improved
performance is primarily associated with correctly finding true lower
multiplicity PVs. Training the same model as that used to produce the red
triangles, but using one of the new KDEs in place of the original KDEs, leads
to worse performance. The algorithm using both of the new KDEs learns how to
combine their information effectively. The plot also shows results for a U-Net
model using the same toy MC as used to train the the “traditional” KDE-to-
hists models. Its performance is better than that of the other models trained
and tested on the same data. Its architecture is discussed below.

##  3 Kernel Density Estimators

Kernel generation converts sparse three-dimensional data into a small number
of feature-rich one-dimensional data sets. Deep neural networks (DNNs) can
transform these into one-dimensional histograms from which PV candidates are
easily extracted. Each of our one-dimensional data sets consists of 4000 bins
along the $z$-direction (beamline), each 100 $\mu$m wide, spanning the active
area of the VELO around the interaction point. In the original KDE, each $z$
bin of the histogram is filled by the maximum kernel value in $x$ and $y$,
where the kernel is defined by

$\mathcal{K}(x,y,z)=\frac{\sum_{\mathrm{tracks}}\mathcal{G}(\mathrm{IP}(x,y)|z)^{2}}{\sum_{\mathrm{tracks}}\mathcal{G}(\mathrm{IP}(x,y)|z)}-\sum_{\mathrm{tracks}}\mathcal{G}(\mathrm{IP}(x,y)|z)\
.$ (1)

In Eqn. (1), $\mathcal{G}(\mathrm{IP}(x,y)|z)$ is a Gaussian function,
centered at $x=y=0$ and evaluated at the impact parameter IP$(x,y)$: the
distance of closest approach of a track projection to a hypothesized vertex at
position $x,y$ for a given $z$. The width/covariance of $\mathcal{G}$ is given
by the IP$(x,y)$ uncertainty/covariance matrix. Finding the maximum
$\mathcal{K}(x,y,z)$ is a two step process. Kernel values are first computed
in a coarse $8\times 8$ grid in $x,y$; the parameters of that search are then
taken as starting points for a MINUIT minimization process to find the maximum
kernel value for each bin of $z$. The values of $x$ and $y$ where the kernel
is maximum are saved as secondary feature sets xMax and yMax

Equation (1) is computed from the location, direction, and covariance matrix
of input tracks. This information can be provided by a standalone toy Monte
Carlo proto-tracking Fang:2019wsd  or from a proposed LHCb Run 3 production
Velo tracking Hennequin:2019itm . While the former uses a heuristic Gaussian
width to estimate the IP$(x,y)$ uncertainty, the latter uses the measured
covariance matrix of the track state closest to beamline to compute the
IP$(x,y)$ covariance. To reduce the computational load, the KDE contributions
are only calculated in regions near points where two tracks pass within
$100\,\mu$m of each other. This is completely safe in terms of finding true PV
positions, but leads to discontinuities in KDE values as a function of $z$.

Figure 3:  This diagram illustrates the deep neural network used to predict an
event’s KDE-A from its tracks’ poca-ellipsoids. Twelve fully connected layers
(all but the last with 50 neurons) populate four 4000-bin channels in the last
of these layers, for each track. These contributions are summed and processed
by two convolutional layers that provide a mechanism for the tracks to
interact. A final fully connected layer was required for the learning to
converge.

Our first attempt to design an algorithm to predict KDEs from track parameters
performed poorly. It was designed to treat each track’s contribution to the
KDE separately, and we guessed that the explicit requirement that a track
contribute to the KDE only when near another track was part of the problem. To
mitigate this, we created a new KDE (denoted KDE-A) based on each track’s
position of closest approach (poca) to the beamline and the error ellipsoid
defined by

$A(\Delta x)^{2}+B(\Delta y)^{2}+C(\Delta z)^{2}+2(D\Delta x\Delta y+E\Delta
x\Delta z+F\Delta y\Delta z)=1$ (2)

where $\Delta x=x-x_{\rm poca}$, etc. At any point in space, a Gaussian
probability associated with this track’s poca-ellipsoid is calculated as

${\cal G}(x,y,z)=\frac{{\rm exp}\left(-0.5(A(\Delta x)^{2}+B(\Delta
y)^{2}+C(\Delta z)^{2}+2(D\Delta x\Delta y+E\Delta x\Delta z+F\Delta y\Delta
z)\right)}{(\mbox{ellipsoid volume}){}^{3/2}}$ (3)

For each voxel, the KDE contribution is calculated as the sum of poca-
ellipsoid probabilities:

${\cal K}_{A}=\sum_{\mathrm{tracks}}{\cal G}(x_{i},y_{i},z_{i})\,.$ (4)

where $i$ denotes the track index. For each bin of $z$, the largest value of
${\cal K}_{A}$ observed in an $(x,y,z)$ voxel is projected out as the KDE-A
value at that value of $z$. However, the performance of the KDE-to-hists
algorithm using KDE-A in place of the original KDE was significantly worse. As
the definition of the original kernel in Eqn. 1 combined linear and quadratic
terms in $\mathcal{G}(\mathrm{IP}(x,y)|z)$, we defined KDE-B using a procedure
parallel to that for defining KDE-A but replacing the linear sum of
probabilities in ${\cal K}_{A}$ with

${\cal K}_{B}=\sum_{\mathrm{tracks}}{\cal G}(x_{i},y_{i},z_{i})^{2}\,.$ (5)

Using both KDE-A and KDE-B as input features, in place of the original KDE,
improved the performance of the KDE-to-hists algorithm, as noted in the
earlier discussion of Fig. 2 (see the red triangles and cyan-filled circles).
This pleasant outcome was not fully anticipated.

As we had hoped, it is easier to teach a DNN to predict KDE-A with reasonably
good fidelity, although achieving satisfactory performance does require a
mechanism that encourages the predictions for each track to interact with the
predictions for the other tracks. We have studied this using toy MC rather
than full LHCb MC as the former data sets are larger than the latter. The
architecture we have used has not been optimized – the choices of the numbers
of layers of each type, the number of channels per layer, the kernel sizes in
the convolutional layers, etc., were chosen almost arbitrarily. With these
caveats, the architecture is illustrated in Fig. 3. Originally, there were 12
fully connected layers that predicted KDE-A contributions for each track
separately and then summed them. The performance of this model was mediocre
(at best). To allow the tracks to interact, the single 4000-bin channel in the
12th layer that was originally meant to be the KDE prediction was replaced by
four latent 4000-bin channels. These are used as input to two convolutional
layers (and one more fully connected layer) that produce a final 4000-bin
prediction. This more complicated model was trained starting with the weights
and biases of the first 11 fully connected layers fixed while those of the new
layers were learned. Once the full model began to display a semblance of
performance, all the model parameters were floated and learned over many
epochs, using progressively larger training samples. Fig. 4 compares KDE-A
(the blue histogram) for the first event in the validation sample with the
corresponding KDE predicted by the model (the red histogram) after about 300
hours of training on an nVidia RTX2080Ti using PyTorch. Qualitatively, the
predicted KDE reproduces the macroscopic features of KDE-A, as seen in the
plots on the left. On a finer scale, there are differences, as seen in the
plot on the right. This shows 50 bins on either side of the location of the
first PV in the event. The predicted KDE peaks at the location of the PV, as
does KDE-A, but its height is lower and its width is greater. This is typical
of the learned KDE, although there is tremendous variation from the region
about one true PV in an event to another and from event to event. A mean-
squared-error function comparing the predicted KDE histogram to the KDE-A is
used to train the model. But the real goal is using a well-trained tracks-to-
KDE model as a step towards defining a tracks-to-hist model that performs at
least as well as the best model in Fig. 2.

Figure 4:  The plot on the left illustrates KDE-A (shown in blue), calculated
as described in the text following Eqn. 5, and the KDE predicted by the DNN
model described in the subsequent paragraph (shown in red) for a single event.
Each of the 4000 bins in the histogram corresponds to a range in the
$z$-direction of $100\,\mu$m. The plot on the right shows a zoomed-in view of
the 50 bins on either side of the location of the second PV in that event,
whose location is denoted by the dashed vertical line.

As a proof-of-principle, we first trained a KDE-to-hists model using KDE-A as
the only input feature set. The efficiency and false positive rate were 94.3%
and 0.27 per event. We then merged the tracks-to-KDE model described above
with this KDE-to-hists model to create a tracks-to-hists model. With the
weights of the combined model set to be those of the independent models, the
efficiency and false positive rate were approximately 88% and 0.88 per event.
We then allowed the combined model to learn improved weights for the tracks-
to-KDE layers while the weights for the KDE-to-hists layers were held constant
and vice versa. In these cases, the improvements were very slow, with the
efficiencies increasing a (small) fraction of a percent and the false positive
rates improving by less than 0.02 over the course of several days training on
nVidia RTX 2080Ti GPUs. Allowing all weights and biases in the combined model
to float produced significantly improved performance. After several days of
training, the efficiency increased to 91.0% and the false positive rate
dropped to 0.80 per event before improvement completely plateaued.

As noted above, the predicted KDEs tend to produce broader features in the
vicinity of true PVs, as seen in Fig. 4. This suggested that we modify the
definitions of efficiency and false positive rate to allow predicted and true
PV positions to be deemed correctly matched when separated by as much as
$750\,\mu$m rather than $500\,\mu$m. With this alternative definition, the
efficiency rises to 93.7% and the false positive rate drops to 0.69 per event.
Although these efficiencies (original and alternative definitions) are not as
high as that of the KDE-to-hists model, and the false positive rates are
higher, the fact that simultaneously learning the combined model parameters
produces significantly better results than simply connecting the two networks
with their independently learned parameters suggests that it should be
possible to create a highly performant tracks-to-hists model. It also suggests
that improving the tracks-to-KDE model so that it produces more narrow peaks
around PV positions will be necessary to significantly improve the merged
model performance.

##  4 Alternative Model Architectures

Figure 5:  This plot illustrates the improved performance of several models
when xMax and yMax are added as perturbative features, as discussed in the
text. For each model, the open markers show the efficiency versus false
positive rate as the cost function asymmetry parameter is varied. The
corresponding filled markers show the performance when the perturbative
features are added.

Figure 6:  This plot illustrates the performance of several models built from
a benchmark model with variations described in the text. The labels for each
model encode the nature of the variations, as described in Table 1 in the
appendix. The number of parameters learned by each model is shown in
parentheses to the right of its label.

The architectures of the models whose performance is discussed in Sec. 2 and
illustrated in Fig. 2 are similar, but vary as discussed at conferences last
year akar2020updated  and above. We have now investigated some of these
variations, and others, more systematically. We have studied a qualitatively
new architecture as well. We find that adding additional the secondary feature
sets xMax and yMax to the models “perturbatively” systematically improves
performance. We also find that models with more parameters generally out-
perform similar models with fewer parameters. From a computer science
perspective, the best performance may be considered optimal. However, when
deploying software for use in a HEP experiment, use of memory and processor
time may be considered as well. For example, it may be appropriate to trade
off small gains in performance if doing so requires only half the computing
resources.

The results of adding perturbative features to three baseline models is
illustrated in Fig. 6. Once a model has been fully trained using only the KDE
feature set, its weights and biases are fixed and a second network is trained
using only the secondary features; the final prediction is the product of the
two networks’ outputs (or generated from a sum or concatenation of their
outputs). After the secondary network is trained in this way, all the weights
and biases of the combined model are trained simultaneously to produce the
final model. We find that initially separating the baseline training from the
secondary feature training is necessary for deep learning to start. In each
case shown (and others not shown to avoid visual congestion), adding xMax and
yMax features to the model perturbatively increases the efficiency and reduces
the false positive rate.

The benchmark model used as the starting point for studies of alternative
models has six convolutional layers, very similar to the model with four
convolutional layers illustrated in Fig. 3 of Ref. Fang:2019wsd . In the
results presented below, xMax and yMax have been added as perturbative
features. The primary goal of these studies was to understand how the fidelity
of a model varies with the number of parameters and the use of skip
connections, batch normalization, etc. Table 1 in the appendix lists the types
of variations that were considered. Results from 5 of these models are shown
in Fig. 6. In general, the larger the number of parameters, the greater the
fidelity of the model. That with the best performance is an adaptation of the
DenseNet architecture huang2018densely  which connects each layer to every
other layer in a feed-forward fashion, the ultimate limit adding skip
connections. Note that these studies were done in parallel with others
reported here. We have not yet incorporated what was learned here into, for
example, the best performing models of Fig. 2.

Figure 7:  This is the architecture of a U-Net model modified for use in our
one-dimensional KDE-to-hists problem. The orange blocks are convolutional
layers followed by ReLU activation functions, the red blocks are down-sampling
operations, and the blue blocks are up-sampling operations. The four arrows
that join down-sampling layers to up-sampling layers of the same
dimensionality denote skip connections.

Inspired by its success solving other tasks that require shape-preservation
between input and output (object detection, semantic segmentation), we have
modified the popular U-Net architecture ronneberger2015unet  for our one-
dimensional problem. A diagram of the architecture is shown in Fig. 7 – orange
blocks are convolutional layers followed by ReLU activation functions, red
blocks are down-sampling operations (max pooling with kernel size and stride
2), and blue blocks are up-sampling operations (transposed convolutions with
kernel size and stride 2). The arrows represent the flow of information
throughout the network. It differs from a traditional autoencoder as some
information bypasses successive layers and is added back in later stages via a
concatenation or addition operation (a skip connection). We investigated this
architecture for a number of reasons. Conceptually, it captures a receptive
field that scales as $2^{n}$, where $n$ is number of pooling operations. This
can be compared to the AllCNN architectures where the receptive field scales
linearly with $n$. We hypothesize that the added ease in down-sampling is
important to achieve the granular resolution we want. In addition, we
hypothesize that the addition of skip connections up to the second-to-last
layer in the network allows the network to “remember” fine-grained detail
about the input, and more effectively localize primary vertices. Experiments
show that removing skip connections one-by-one consistently degrades results
by a minor amount. Including all four skip connections, the architecture
achieves an efficiency of 95% with a false-positive rate of 0.16 per event,
slightly better than the best AllCNN architecture using the same toy MC data
and features. We anticipate that this architecture will perform at least as
well as the best AllCNN architecture when trained using full LHCb MC plus
KDE-A and KDE-B as input features.

##  5 Summary and Conclusions

Since we presented results at Connecting-the-Dots 2020, we have tested AllCNN
KDE-to-hists models using full LHCb Monte Carlo rather than toy Monte Carlo
and find they are (i) more performant without re-training and (ii) even more
performant with re-training. If we replace the original KDE with a pair of
KDEs, each built from individual track contributions with no pairwise
interactions used explicitly, the performance is better again. We have built a
deep learning tracks-to-KDE algorithm that reproduces the new KDE of Eqn. 5
with sufficient fidelity that a merged tracks-to-hist algorithm achieves
almost as high efficiency as the underlying KDE-to-hist algorithm, albeit with
worse resolution and a significantly higher false positive rate. We conclude
that the overall approach is sound, but we need to construct and train a
better tracks-to-KDE algorithm.

We have also performed systematic studies of alternative AllCNN models and a
very different model inspired by U-Net. The results of these studies indicate
that the capacity of a model to learn increases with the number of parameters
and they suggest specific approaches for improving the performance of the best
KDE-to-hist model so far. The very similar performance of the model inspired
by U-Net and the best AllCNN model trained on the same KDE invites the
question of how two such different architectures have “learned” the same
“concepts” and whether such “understanding” can be generalized.

Two physics insights/hypotheses underly our current approach to developing
hybrid deep learning algorithms for identifying and locating primary vertices.
First, tracks from primary vertices pass close enough to the beamline that
each can be characterized by its poca-ellipsoid. Second, the sparse three-
dimensional information encoded in these ellipsoids can be transformed to
rich, one-dimensional KDEs that can, in turn, serve as input feature sets for
relatively simple deep neural networks. These can be trained on large samples
of simulated events (hundreds of thousands) to produce one-dimensional
histograms that can be processed by traditional, heuristic algorithms to
extract information for the next stage of event reconstruction and
classification.

Our specific algorithms are being designed for use in LHCb with its Run 3
detector. As the salient characteristics exist in the ATLAS and CMS
experiments, as well, appropriately modified versions might work for them. A
heuristic ATLAS vertexing algorithm ATLAS:2019jmx ; Sanderswood:2019rta ,
already uses some conceptually similar ideas, as does an algorithm
Bocci:2020pmi  designed for the CMS detector, upgraded for operation in the
high luminosity LHC era. As usual, the devil is in the details. Developing and
deploying machine learning inference engines that are highly performant and
satisfy computing system constraints requires sustained effort. The results
reported here should encourage work focussed on using deep neural networks for
identifying vertices in high energy physics experiments.

##  6 Acknowledgments

The authors thank the LHCb computing and simulation teams for their support
and for producing the simulated LHCb samples used in this paper. The authors
also thank the full LHCb Real Time Analysis team, especially the developers of
the VELO tracking algorithm Hennequin:2019itm  used to generate the “full LHCb
MC” results presented in Fig 2.

This work was supported, by the U.S. National Science Foundation under
Cooperative Agreement OAC-1836650 and awards PHY-1806260, OAC-1739772, and
OAC-1740102.

## References

  * (1) R. Aaij et al. (LHCb), _LHCb Trigger and Online Upgrade Technical Design Report_ (2014), https://cds.cern.ch/record/1701361
  * (2) F. Reiss, _Excerpts from the LHCb cookbook â RECEPTS for testing lepton universality and reconstructing primary vertices_ (2020), presented 18 12 2020, https://cds.cern.ch/record/2749592
  * (3) F. Reiss et al. (LHCb), _Fast parallel Primary Vertex reconstruction for the LHCb Upgrade_ (2020), talk presented at Connecting-the-Dots 2020, https://indico.cern.ch/event/831165/timetable/?view=standard
  * (4) R. Aaij et al. (LHCb), JINST 14, P04013 (2019), arXiv:1812.10790
  * (5) R. Fang, H.F. Schreiner, M.D. Sokoloff, C. Weisser, M. Williams, J. Phys. Conf. Ser. 1525, 012079 (2020), arXiv:1906.08306
  * (6) S. Akar, T.J. Boettcher, S. Carl, H.F. Schreiner, M.D. Sokoloff, M. Stahl, C. Weisser, M. Williams, _An updated hybrid deep learning algorithm for identifying and locating primary vertices_ (2020), arXiv:2007.01023
  * (7) A. Hennequin, B. Couturier, V. Gligorov, S. Ponce, R. Quagliani, L. Lacassagne, JINST 15, P06018 (2020), arXiv:1912.09901
  * (8) O. Ronneberger, P. Fischer, T. Brox, _U-Net: Convolutional Networks for Biomedical Image Segmentation_ , in _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_ , edited by N. Navab, J. Hornegger, W.M. Wells, A.F. Frangi (Springer International Publishing, Cham, 2015), pp. 234–241, ISBN 978-3-319-24574-4, arXiv:1505.04597
  * (9) G. Huang, Z. Liu, L. Van Der Maaten, K.Q. Weinberger, _Densely Connected Convolutional Networks_ , in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (2017), pp. 2261–2269, arXiv:1608.06993
  * (10) _Development of ATLAS Primary Vertex Reconstruction for LHC Run 3_ (2019), ATL-PHYS-PUB-2019-015, https://cds.cern.ch/record/2670380
  * (11) I. Sanderswood (ATLAS), _Development of ATLAS Primary Vertex Reconstruction for LHC Run 3_ , in _Connecting the Dots and Workshop on Intelligent Trackers_ (2019), arXiv:1910.08405
  * (12) A. Bocci, M. Kortelainen, V. Innocente, F. Pantaleo, M. Rovere, Front. Big Data 3, 601728 (2020), arXiv:2008.13461

## Appendix: Models and Tags

The models discussed in the studies of systematic variations are derived from
a benchmark AllCNN model. The relationships between these models are shown in
Fig. 8. As indicated in Table 1, the numbers of layers, channels per layer,
use of skip connections, batch normalization, etc., could be varied. Not all
combinations were studied, but the conclusion presented in the body of the
text – that increasing the numbers of parameters increases the capacity of a
network to learn well, almost independently of their origins – seems to be
robust.

Figure 8:  This figure shows how models studied relate to each other. The tags
indicating the variations in the architectures are described in Table 1.

These are the tags used to label each model name. The tags determine some of
the key features of the model at a high level.

Table 1: Neural network model tags for names Tag | V (Meaning)
---|---
#S | number of skip connections, if any
#L | number of layers
BN | Batch Normalization
ACN | AllCNN “family” of models
#i# | step number out of total steps
RC# | reduced channel size, followed by iteration number
IC# | increased channel size, followed by iteration number
RK# | reduced kernel size, followed by iteration number
IK# | increased kernel size, followed by iteration number
C | concatenation of perturbative and non-perturbative layers at end
BM | benchmark; changes made to future models based on tagged reference model

The tag hierarchy is of the format:

BM-ACN-#i#-P-#L-#S-BN-RC-#-IC#-RK#-IK#-C