Vacuum-field-induced state mixing

Diego Fernández de la Pradilla1$\star$,
Esteban Moreno1 and
Johannes Feist1

1 Departamento de Física Teórica de la Materia Condensada and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain

${}^\star$<EMAIL_ADDRESS>
August 27, 2024

§ ABSTRACT

By engineering the electromagnetic vacuum field, the induced Casimir-Polder shift (also known as Lamb shift) and spontaneous emission rates of individual atomic levels can be controlled. When the strength of these effects becomes comparable to the energy difference between two previously uncoupled atomic states, an environment-induced interaction between these states appears after tracing over the environment. This interaction has been previously studied for degenerate levels and simple geometries involving infinite, perfectly conducting half-spaces or free space. Here, we generalize these studies by developing a convenient description that permits the analysis of these non-diagonal perturbations to the atomic Hamiltonian in terms of an accurate non-Hermitian Hamiltonian. Applying this theory to a hydrogen atom close to a dielectric nanoparticle, we show strong vacuum-field-induced state mixing that leads to drastic modifications in both the energies and decay rates compared to conventional diagonal perturbation theory. In particular, contrary to the expected Purcell enhancement, we find a surprising decrease of decay rates within a considerable range of atom-nanoparticle separations. Furthermore, we quantify the large degree of mixing of the unperturbed eigenstates due to the non-diagonal perturbation. Our work opens new quantum state manipulation possibilities in emitters with closely spaced energy levels.

§ INTRODUCTION

It is well known that atomic properties are modified due to the interaction with the quantized electromagnetic (EM) vacuum field supported by macroscopic bodies <cit.>. In the weak-coupling regime, this changes both the atomic linewidths (Purcell effect) <cit.> and energies (Lamb or Casimir-Polder [CP] shifts) <cit.>. These modifications have wide-ranging applications in fields such as optics or atomic and soft matter physics, including the design of efficient single photon sources <cit.>, the atomic force microscope <cit.>, new atom trapping methods <cit.> or the precise manipulation of atomic properties with tunable nanostructures <cit.>. Theoretical descriptions of these effects are commonly perturbative, using either standard perturbation theory or open quantum systems approaches <cit.>, although efforts to go beyond the purely perturbative regime have also been published <cit.>. When the interactions are weak, the effect of the environment is customarily treated for each atomic state independently, giving rise to simple diagonal energy shifts and decay rates. However, for subsets of near-degenerate atomic states, the CP shift and/or spontaneous emission rates may be of the same order as the energy differences within the subset, suggesting that the above treatment is not consistent, even if the light-matter coupling is perturbative. This has been discussed in the literature for atoms in free space <cit.>.

In this work, we show that the standard diagonal perturbation approach indeed fails when field-induced shifts are comparable to the energy level differences, requiring the treatment of environment-induced interactions between the levels <cit.>. Recently, this issue has been tackled with an open quantum systems' framework designed for structures with closely spaced levels <cit.>. In that work, the standard Bloch-Redfield equation <cit.> was turned into a Lindblad equation <cit.>, with the corresponding benefit of guaranteed positive populations, while simultaneously eluding the usual secular approximation that neglects the couplings between non-degenerate states <cit.>. Here, we extend this framework to incorporate the effect of the counter-rotating terms in the light-matter Hamiltonian and construct a master equation that accurately represents the off-diagonal CP and decay terms, which we expect to be relevant in any system with subsets of near-degenerate levels. From the Lindblad equation, we extract an effective non-Hermitian Hamiltonian that determines the dynamics of a subset of levels and in turn enables a quantitative exploration of the vacuum-field-induced state mixing. We illustrate the effects of the off-diagonal terms by applying the above steps to a system comprised of a hydrogen atom close to an aluminum nitride (AlN) nanoparticle (NP), and study the impact of the off-diagonal couplings on the dynamics of the atom. We find strong modifications to the level structure and observe significant state mixing at atom-NP separations on the order of \(100~\)nm. Consequently, the atomic dynamics close to the NP cannot be understood without consideration of the effects discussed in this work. The situation treated here lies between conventional weak coupling (where light-matter interactions can be treated perturbatively and states can be considered independently) and strong coupling (where light and matter excitations mix significantly due to non-perturbative interactions). In this novel regime of “strong weak coupling”, perturbative light-matter interactions lead to significant state mixing within the matter component.

§ METHODS

§.§ Macroscopic quantum electrodynamics

We describe the interaction between atoms and the EM field supported by macroscopic bodies within macroscopic quantum electrodynamics (MQED) <cit.>. The corresponding Power-Zienau-Woolley light-matter Hamiltonian in the dipole approximation <cit.> is

\begin{equation} \label{eq.Hamiltonian}
H = H_\textrm{at} + H_\textrm{f} - \textbf{d}\cdot\textbf{E}(\mathbf{r}_\mathrm{at}).
\end{equation}

Here, \(H_\textrm{at}\) is the matter Hamiltonian, emphasizing that we treat a single atom. The field Hamiltonian,

\begin{equation}
H_\mathrm{f} = \sum_\lambda\int\mathrm{d}^3r\int\textrm{d}\omega\,\hbar\omega\,\mathbf{f}_\lambda^\dagger(\mathbf{r},\omega)\cdot \mathbf{f}_\lambda(\mathbf{r},\omega),
\end{equation}

contains the (bosonic) polaritonic annihilation and creation operators \(\mathbf{f}_\lambda(\mathbf{r},\omega)\) and \(\mathbf{f}_\lambda^\dagger(\mathbf{r},\omega)\) that describe both the purely electromagnetic and the macroscopic polarization fields. Here, the index \(\lambda = \{e, m\}\) labels the electric or magnetic nature of the excitations, and the integrals are over all space and over all positive frequencies. The last term is the dipolar interaction between the atom with dipole operator \(\mathbf{d}\) and the electric field \(\mathbf{E}(\mathbf{r})\) evaluated at the atomic position \(\mathbf{r}_\mathrm{at}\), where

\begin{equation}
\mathbf{E}(\mathbf{r}) = \sum_\lambda\int\mathrm{d}^3s\int \mathrm{d}\omega \, \mathbf{G}_\lambda(\mathbf{r},\mathbf{s},\omega)\cdot\mathbf{f}_\lambda(\mathbf{s},\omega) + \mathrm{H.c.},
\end{equation}

\begin{align*}
\mathbf{G}_e(\mathbf{r},\mathbf{s},\omega) &= i \frac{\omega^2}{c^2} \sqrt{\frac{\hbar}{\pi\varepsilon_0}\Im{\varepsilon(\mathbf{s},\omega)}}\mathbf{G}(\mathbf{r},\mathbf{s},\omega), \\
\mathbf{G}_m(\mathbf{r},\mathbf{s},\omega) &= i\frac{\omega}{c}\sqrt{\frac{\hbar}{\pi\varepsilon_0}\frac{\Im{\mu(\mathbf{s},\omega)}}{|\mu(\mathbf{s},\omega)|^2}} \Big[\nabla_\mathbf{s}\times\mathbf{G}(\mathbf{s},\mathbf{r},\omega)\Big]^T.
\end{align*}

Here, \(\varepsilon\) and \(\mu\) stand for the electric and magnetic response functions, respectively. \(\mathbf{G} = \mathbf{G}^0+\mathbf{G}^\mathrm{scatt}\) is the classical electromagnetic Green tensor, separated in its free-space and scattering contributions. In the weak-coupling regime, \(\mathbf{G}^0\) is responsible for the free-space Lamb shift <cit.>, a NP-independent contribution that can be simply reabsorbed in \(H_\mathrm{at}\). Compared to the effects we study in this work, this is a negligible correction that we discard entirely in the following.

(a) Sketch of the system. (b) Simplified level structure of the hydrogen atom with Bohr levels and fine structure splitting (not to scale). Diagonal lines: dipolar transitions allowed from the \(n=4\) Bohr level to the \(n=5\) and \(n=3\) Bohr levels. (c) Top: spectral density of the AlN NP, obtained by setting \(d=50\) nm in <ref>, bottom: result of the integral in <ref>. \(e\) is the absolute value of the charge of the electron and \(a_0\) is the Bohr radius. Vertical line: transition frequencies of the atom from \(n=7\) to \(n=6\).

For concreteness, we focus on a hydrogen atom interacting with a spheroidal AlN NP (see <ref>). It should be noted, however, that the following arguments are of broader generality and applicable to a wide range of physical systems, provided that the energies and transition dipole moments of the atom and the Green tensor of the nanostructure are accessible, a general requirement of CP calculations. We choose this NP shape and material for two reasons: (i) the EM resonances along the symmetry axis (\(z\)) enhance the atomic transitions mediated by \(E_z\) with respect to the other components, and (ii) the energy range of the EM resonances coincides with the hydrogenic transition we want to target. Hence, this system provides a realistic and not overly complicated testing ground for the formalism developed, and will allow us to illustrate the effects of the off-diagonal vacuum shifts.

In <ref>, \(H_\mathrm{at}\) is diagonal, and its eigenvalues include fine structure corrections <cit.>:

\begin{equation}
E_{nj} = \left(-\frac{1}{2n^2} - \frac{\alpha^2}{2n^3} \left[ \frac{1}{j+\frac{1}{2}} - \frac{3}{4n}\right] \right)E_\textrm{h},
\end{equation}

where \(n\), \(j\), \(\alpha\) and \(E_\mathrm{h} \simeq 27.2~\)eV are the main quantum number, the total electronic angular momentum, the fine structure constant and the Hartree energy, respectively. The energy levels, schematically shown in <ref>, are distributed in well-separated Bohr levels labeled by \(n\), corrected with the fine structure splitting \(\Delta_F\), a \(j\)-dependent quantity that is 4 or more orders of magnitude smaller. This energy scale is small enough that the CP induced interaction between fine structure states with the same \(n\) can be relevant. The NP is a spheroid with major and minor axes \(M=200\) nm and \(m=140\) nm, with the AlN dielectric permittivity taken from <cit.>. For this NP, the phonon-polariton resonances lie close to the transition energy between the hydrogenic \(n=7\) and \(n=6\) states. Specifically, we will focus on the off-diagonal effects within the \(n=7\) level for an atom located along the symmetry axis \(z\) of the NP.

§.§ Master equation and effective non-Hermitian Hamiltonian

To describe the dynamics of the field-modified atomic levels and their mixing, we derive a Lindblad equation for the atomic density matrix \(\rho\) by considering the EM fields as a weakly coupled bath and perturbatively tracing out the EM degrees of freedom. We start from the standard open quantum systems approach, which leads to the so-called Bloch-Redfield equation <cit.>:

\begin{align} \nonumber
\dot{\rho} = - \frac{i}{\hbar}\left[H_\mathrm{at},\rho\right]+ \sum_{abcd}\Bigg[ &- i\Big( \Lambda_{ca,db}(\omega_{bd}) \ket{a}\braket{c}{d}\bra{b}\rho - \Lambda_{ca,db}(\omega_{ac}) \rho \ket{a}\braket{c}{d}\bra{b}  \Big) \\ \nonumber
&+ i\big[\Lambda_{ca,db}(\omega_{bd})- \Lambda_{ca,db}(\omega_{ac}) \big] \ketbra{d}{b}\rho\ketbra{a}{c}   \\ \nonumber
&- \frac{1}{2}\Big( \Gamma_{ca,db}(\omega_{bd}) \ket{a}\braket{c}{d}\bra{b}\rho + \Gamma_{ca,db}(\omega_{ac}) \rho \ket{a}\braket{c}{d}\bra{b} \Big)\\ \label{eq.blochredfield_rot_inv}
&+ \frac{1}{2}\big[\Gamma_{ca,db}(\omega_{bd}) + \Gamma_{ca,db}(\omega_{ac}) \big]\ketbra{d}{b}\rho\ketbra{a}{c} \Bigg],
\end{align}

where \(\rho\) is the atomic density matrix, the Latin indices $a,b,c,d$ denote atomic eigenstates and \(\omega_{ab} = (E_a-E_b)/\hbar\). The rotationally invariant quantities \(\Gamma_{ca,db}(\omega)\) and \(\Lambda_{ca,db}(\omega)\) are given by

\begin{align}
\Gamma_{ca,db}(\omega) &= \mathbf{d}_{ca}^*\cdot \boldsymbol{\gamma}(\omega) \cdot  \mathbf{d}_{db} \\
\Lambda_{ca,db}(\omega) &= \mathbf{d}_{ca}^*\cdot \boldsymbol{\lambda}(\omega) \cdot  \mathbf{d}_{db},
\end{align}

where \(\mathbf{d}_{ca}\) is a matrix element of the atomic dipole operator, and \(\boldsymbol{\gamma}\) and \(\boldsymbol{\lambda}\) are defined as

\begin{align}
\boldsymbol{\gamma}(\omega) &= 2\pi \mathbf{J}(\omega),  \label{eq.gammaab} \\
\boldsymbol{\lambda}(\omega) &= \mathcal{P}\int \mathrm{d}\omega' \frac{\mathbf{J}^\mathrm{scatt}(\omega')}{\omega-\omega'}.\label{eq.lambdaab}
\end{align}

Here, \(\mathcal{P}\) denotes the principal value and \(\mathbf{J}(\omega)\) is the spectral density of the EM field,

\begin{equation} \label{eq.spectral_density}
\mathbf{J}(\omega) = \frac{\omega^2}{\hbar\pi\epsilon_0c^2} \Im{\mathbf{G}(\mathbf{r}_\mathrm{at},\mathbf{r}_\mathrm{at},\omega)}.
\end{equation}

We have replaced \(\mathbf{G}\) by \(\mathbf{G}^\mathrm{scatt}\) in \(\mathbf{J}^\mathrm{scatt}\) (<ref>), since the free-space contribution is assumed to be included in the bare atomic Hamiltonian, and it is also smaller than the NP-induced effects discussed in this work. We note that while we indicate the complex conjugation of the dipole matrix elements, it is possible to choose an atomic basis in which they are real; thus, both \(\Gamma\) and \(\Lambda\) would be real quantities. Note that the expressions for \(\boldsymbol{\lambda}\) contain both the so-called resonant contributions, which can be shown to be proportional to \(\Re{\mathbf{G}^\mathrm{scatt}}\), and non-resonant contributions to the energy shift <cit.>. The electromagnetic Green tensor is computed using the boundary element method implemented in SCUFF-EM <cit.>.

<ref> does indeed describe the bath-induced interaction between levels, but allows for non-physical dynamics, as it does not guarantee the positivity of the density matrix. A standard secularization procedure leads to a completely positive Lindblad equation, but removes the crucial off-diagonal terms describing state mixing (for more details see <ref> of the appendix). Instead of secularization, we extend the approach of Ref. <cit.> to obtain a completely positive Lindblad equation for near-degenerate levels by including the effect of the counter-rotating terms of the dipolar interaction. This approach consists in replacing both \(\Lambda_{ca,db}(\omega_{bd})\) and \(\Lambda_{ca,db}(\omega_{ac})\) with their geometric mean \(\tilde{\Lambda}_{ca,db} = \sqrt{\Lambda_{ca,db}(\omega_{bd})}\sqrt{\Lambda_{ca,db}(\omega_{ac})}\), and the same for \(\Gamma_{ca,db}(\omega_{bd})\) and \(\Gamma_{ca,db}(\omega_{ac})\). Similar ideas have also been proposed elsewhere in the literature <cit.>. When applied to <ref>, this replacement symmetrizes the pairs of indices \(ca\) and \(db\). Then, for the symmetric geometry of <ref> where \(\mathbf{G}\) is a diagonal matrix, the resulting master equation can be rewritten as

\begin{align}\label{eq.master}
\dot{\rho} = \frac{-i}{\hbar} \left[H_\mathrm{at} + H_\mathrm{CP},\rho\right] + \sum_{\delta, n} L_{\Sigma_\delta^{(n)}}[\rho].
\end{align}

Here, \(H_\mathrm{CP}\) is the CP shift, \(\Sigma_\delta^{(n)}\) are decay operators, and \(L_A[\rho] = A\rho A^\dagger -\frac{1}{2}\{A^\dagger A,\rho\}\) is a Lindblad dissipator. There is a decay operator for each spatial component in the spherical basis \(\delta = 0,\pm1\), and for each Bohr level \(n\). The CP shift is given by \(H_\mathrm{CP}=\hbar\sum_{\delta, n} D^{(n)\dagger}_\delta D^{(n)}_\delta\), where

\begin{align}
\bra{j}D^{(n)}_\delta\ket{n} &=\sqrt{\lambda_{\delta\delta}(\omega_{nj})} d^\delta_{jn},
\end{align}

and the decay operators can be expressed as

\begin{align}
\bra{j}\Sigma_\delta^{(n)}\ket{n} &= \sqrt{\gamma_{\delta\delta}(\omega_{nj})} d^\delta_{jn}.
\end{align}

In the definitions above, the atomic states \(\ket{n}\) and \(\ket{j}\) belong to the \(n\)th and \(j\)th Bohr level, respectively. Had the counter-rotating terms in the dipolar coupling not been taken into account, the Hermitian dipole operators in <ref> would have been replaced by raising and lowering operators, such that only terms where \(\omega_j \leq \omega_n\) would be present in <ref>. Instead, our description also incorporates the CP shift contribution given by states with \(\omega_j > \omega_n\). More details about the derivation steps can be found in <ref> of the appendix.

There are two details that require addressing to completely justify the above Lindblad master equation. The first one relates to the free-space Lamb shift. Although we have done the previous manipulations already assuming that \(\mathbf{G}^\mathrm{scatt}\) is the only relevant part of the energy shift (see <ref>), strictly speaking there is also the free-space part. Even if the free-space Lamb shift is indeed negligible compared to the NP-assisted shift, the geometric mean introduces a cross-term between free-space and scattered contributions that one could expect to be larger than the pure free-space shift. However, if \(\lambda_i = \lambda_{i,0} + \lambda_{i,s} \) then the shift can be expanded as \(\sqrt{\lambda_1 \lambda_2} \approx \bar{\lambda}_{s} + \frac12 \left(\frac{\bar{\lambda}_{s}}{\lambda_{1,s}} \lambda_{1,0} + \frac{\bar{\lambda}_{s}}{\lambda_{2,s}} \lambda_{2,0} \right)\), where \(\bar{\lambda}_{s} = \sqrt{\lambda_{1,s} \lambda_{2,s}}\). Since the dynamically relevant terms are those where \(\lambda_{1,s} \approx \lambda_{2,s} \approx \bar{\lambda}_{s}\), the final correction due to the free-space contribution turns out to be \(\approx \frac12(\lambda_{1,0}+\lambda_{2,0})\), that is, of the same order as the free-space shift itself. Therefore, the simple argument to discard the free-space contribution turns out to be enough for the cross-terms that appear with the geometric mean too. Finally, we note that for an atom interacting with a thermal bath, the expression in <ref> would contain additional terms proportional to the occupation number \(n_T(\omega) = (e^{\hbar\omega/k_BT} - 1)^{-1}\), but at the relevant transition frequencies this can be shown to be very small for laboratory temperatures. Therefore, these thermal contributions have been neglected.

While solving <ref> is considerably more affordable than a direct solution of the Schrödinger equation with the Hamiltonian from <ref>, several faithful approximations allow for further simplification, succinctly described below, with more details and an explicit check of the validity of the approximations given in <ref> and <ref>, respectively. First, for the dynamics within a single Bohr level, here \(n=7\), we can discard the states with \(n\neq7\) and write a closed set of equations for \(n=7\), due to the large difference in the energy scales associated to the Bohr transition energies and the environment-induced perturbations. The other states are then considered only implicitly as intermediate virtual states that contribute to the CP and decay terms. Furthermore, within this subspace, the “quantum jump” terms \(\Sigma_\delta^{(7)} \rho \Sigma_\delta^{(7)\dagger}\) in <ref> are negligible since they are proportional to \(J_{\delta\delta}(\Delta_F)\), and the spectral density approaches zero for small frequencies. With these approximations, the dynamics within the \(n=7\) subspace is described by an effective non-Hermitian Hamiltonian

\begin{equation} \label{eq.effHamiltonian}
H_\mathrm{eff}^{(7)} = H_\mathrm{at}^{(7)} + \hbar\sum_\delta \left( D_\delta^{(7)\dagger} D_\delta^{(7)} - \frac{i}{2} \Sigma_\delta^{(7)\dagger}\Sigma_\delta^{(7)} \right),
\end{equation}

where \(H_\mathrm{at}\) has been projected onto the \(n=7\) subspace. Last, due to the axial symmetry of the system (see <ref>), the \(z\) component of the total angular momentum is conserved and <ref> consists of independent blocks for each value of \(m_j\). In the following, we focus on the subspace \(m_j=1/2\), which reduces the number of states to be considered to \(7\) for this particular case.

§ RESULTS

The above derivation significantly simplifies the analysis of the dynamics. In particular, the effective Hamiltonian can be diagonalized, and the real and imaginary parts of its eigenvalues correspond to the energies and decay rates of the states including the vacuum-field-induced state mixing. These energies and decay rates are shown in <ref> and <ref>, respectively. Since the CP shift is dominated by an overall attraction to the surface (inset of <ref>), we plot it relative to the average value for each separation \(d\), revealing a completely different and much more complex structure compared to the fully secularized, diagonal model. In particular, clear avoided crossings highlight the relevance of the off-diagonal terms, similar to effects found in interatomic interactions in <cit.>, but here occurring within a single atom. Similarly striking differences between both models appear in the decay rates shown in <ref>. Due to the off-diagonal terms, the decay rates cross each other several times. Particularly prominent is the vacuum-field-induced generation of a state that becomes more protected against spontaneous decay as the atom approaches the NP for separations between about \(d = 40\) nm and \(60\) nm. This is in stark contrast to the behavior when the states are treated independently, for which the effect of quenching leads to monotonic increase of the Purcell factor and thus decay rate with decreasing separation <cit.>. Here, the same environment that induces the decay also produces the interactions that mix the states and leads to the formation of a state protected from the influence of the environment. We note that the subradiant or metastable state created by field-induced mixing has a smaller decay rate than any of the original eigenstates of the atom at the same distance when mixing is not included.

Every quantity is plotted against the atom-NP separation \(d\). (a) Energies of the eigenstates of the effective non-Hermitian Hamiltonian. (b) Participation ratios \(P_i\) indicating the degree of mixing of the eigenstates. (c) Decay rates of the eigenstates of the effective non-Hermitian Hamiltonian. (d) Decay rates with the NP shape chosen to optimize the decay rate reduction: \(M=200\) nm and \(m=120\) nm. For plots (a), (c) and (d), the solid black lines correspond to the full model with off-diagonal terms, and the green dashed lines to the model without the off-diagonal terms. For plot (b), the green lines are the participation ratio of each eigenstate, and the black line represents the average participation ratio.

The emergence of this protected state can be understood by realizing that the system approaches an idealized situation in which one of the states in the \(n=7\), \(m_j = 1/2\) manifold is fully decoupled from the EM environment. This situation would occur if we could ignore the (i) fine structure, (ii) \(x\)- and \(y\)-polarized electric fields in <ref>, and (iii) contributions of states outside of the \(n=6\) Bohr level to the CP shift and spontaneous decay. Then, \(7\) states in \(n=7\) couple to \(6\) states in \(n=6\) through the single operator \(d_z\), and there is always one superposition (“dark state”) with vanishing coupling. For the realistic system, this idealized situation is approached for various reasons. First, the elongated shape of the NP suppresses the $xx$ and $yy$ components of $J$ and $\lambda$ compared to the $zz$ component. Second, the coincidence of the first peak of $J_{zz}$ and $\lambda_{zz}$ with the energy of the transition from $n=7$ to $n=6$ enhances the contributions from \(n=6\) intermediate states compared to other Bohr levels. Lastly, when the CP shifts become greater than the fine structure, the latter becomes a perturbative correction that can be neglected to lowest order. Based on these considerations, we change the aspect ratio of the NP by decreasing the minor axis to \(m=120\) nm in order to amplify the protection of the state. As shown in <ref>, the minimum decay rate becomes an order of magnitude smaller than the naive expectation without off-diagonal terms, unambiguously demonstrating that the off-diagonal terms can significantly impact the structure of the atom and cannot be neglected in a realistic description.

Next, we evaluate the amount of vacuum-field-induced state mixing. The eigenstates \(\ket{\psi}\) of <ref> are linear superpositions of the fine structure basis states \(\ket{\phi_k}\), and the degree of this mixing can be quantified using the so-called participation ratio \(P\) <cit.>, defined as

\begin{equation} \label{eq.participation}
P(\ket{\psi}) = \left[\sum_k |\braket{\psi}{\phi_k}|^4 \right]^{-1}.
\end{equation}

It measures the number of basis states “equally” contributing to the normalized state \(\ket{\psi}\), with possible values ranging from 1 to the number of basis states (7 for the case studied here). For example, for a state of the form \(\ket{\psi} = \sqrt{1/n}\sum_{k=1}^n e^{i\theta_k}\ket{\phi_k}\), \(P\) equals \(n\). In <ref>, we show the participation ratio of the eigenstates of <ref> as a function of the atom-NP separation \(d\), with green lines showing \(P\) for each eigenstate and the thick black line representing the average over all states. We find state mixing to be negligible for separations above \(\approx150\) nm, indicating that the off-diagonal contributions to <ref> are too small to effectively couple the states. For shorter distances, the magnitude of the off-diagonal terms, \(\big|\braOket{i}{H_\mathrm{eff}}{j}\big|\), becomes comparable to the difference of the corresponding diagonal elements, \(\big|\braOket{i}{H_\mathrm{eff}}{i} - \braOket{j}{H_\mathrm{eff}}{j}\big|\), and the states mix appreciably. In particular, clear peaks in \(P\) appear when the diagonal CP shifts bring initially detuned states into resonance, such that the off-diagonal elements dominate more easily. Despite the non-monotonic behavior, overall the participation ratios tend to increase as the atom approaches the NP until they stabilize at about \(d=30\) nm. At closer distances, vacuum-field couplings determine the eigenstates and dominate over the fine structure, which becomes a small perturbation of these new eigenstates. We note that for the setup studied here, almost complete mixing of some atomic states is achieved, with values of \(P\) larger than 5, close to the theoretical maximum of 7.

§ CONCLUSION

In conclusion, we have shown that vacuum-field-induced interactions can significantly mix groups of near-degenerate levels in atoms and must, therefore, be included to accurately characterize the dynamics. We have derived a completely positive Lindblad master equation that provides a precise description of these situations and allows an interpretation of the field-modified atomic structure in terms of an effective non-Hermitian Hamiltonian. For concreteness, we have applied these general ideas to a hydrogen atom coupled to an AlN NP. This leads to striking new features in the atomic structure, such as avoided level crossings and, surprisingly, a decrease of the decay rate for a particular eigenstate with decreasing distance to the NP even though the Purcell factor for each uncoupled state grows monotonically. This illustrates that the off-diagonal terms can even have counterintuitive consequences. Deeper exploration of the eigenstates reveals that the atomic structure in this regime greatly differs from the original fine structure of the free-space atom, even though the atom-field interactions are perturbative and the atom remains a well-defined entity, in contrast to the strong light-matter coupling regime. From an atomic physics' perspective, the hydrogen atom treated here becomes “unrecognizable” as the atomic structure and spectroscopic properties within each sublevel change completely. We note that while we treat a specific setup here, the framework can be straightforwardly applied to other nanostructures (e.g., graphene <cit.>) and emitters or level splittings (e.g., due to hyperfine structure in Rydberg atoms <cit.>). Our work thus extends the regime where vacuum-field-induced forces and decay rates are accurately described and opens the door to new strategies for developing quantum state manipulation platforms based on off-diagonal vacuum-induced effects.

§ ACKNOWLEDGEMENTS

Funding information
This work has been funded by the Spanish Ministry of Science, Innovation and Universities-Agencia Estatal de Investigación through the FPI contract No. PRE2019-090589 as well as grants RTI2018-099737-BI00, PID2021-125894NB-I00, and CEX2018-000805-M (through the María de Maeztu program for Units of Excellence in R&D). We also acknowledge financial support from the Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM) of the Community of Madrid, and from the European Research Council through grant ERC-2016-StG-714870.

§ DERIVATION OF THE MASTER EQUATION

The derivation of the Lindblad master equation used in this work, <ref>, is closely based on the one presented in <cit.>, modified to include the counter-rotating terms of the light-matter Hamiltonian and a realistic electromagnetic environment, with three spatial components and non-trivial structure. We revisit the derivation here and highlight the additions and differences compared to <cit.>. For simplicity, the derivation is presented in a way that directly relates to the illustrative physical system of the main text, that is, a hydrogen atom. However, this is not a true limitation of the approach and approximations, as long as one considers level structures with distinct subsets of closely spaced states, a common feature of atomic systems due to fine structure or hyperfine structure splittings. We start by describing features of the Bloch-Redfield (BR) equation, and then explain the customary secular approximation, a procedure known to yield a completely positive Lindblad master equation. This equation would systematically neglect the off-diagonal terms discussed in this work. Then, we take the BR equation and perform a series of approximations that lead to a different Lindblad equation, maintaining the non-secular terms.

§.§ Comments on the Bloch-Redfield master equation

The BR equation for our system is given by <ref>. Although it is rather complicated, the physical interpretation of each line is simple: the first two lines of the sum are responsible for the Casimir-Polder (CP) shifts, while the third and fourth lines describe decay processes. While the CP terms can often be neglected, this is not the case for the system we study, as they are of the same order as the hydrogenic fine structure. For the same reason, we must include the counter-rotating (CR) terms in the full light-matter Hamiltonian. Otherwise, \(\boldsymbol{\lambda}(\omega)\) would only be evaluated at non-negative frequencies and <ref> would miss significant contributions to the CP terms arising from the negative frequencies. The CR terms only affect the energy shift, as the decay terms \(\boldsymbol{\gamma}(\omega)\) vanish at negative frequencies. It is worth noting that even without considering the CR terms, the equation already includes the basis for the off-diagonal CP terms we discuss in the main text, albeit in a complex manner that is hard to disentangle.

The BR equation has several drawbacks: First, it does not guarantee positivity of the density matrix. Although it has been shown that these deviations from physical density matrices are negligible when the approximations made in deriving the BR equation are valid <cit.>, dealing with formally unphysical density matrices requires additional care. Second, the BR matrix is characterized by a superoperator of dimension $N^2 \times N^2$, where $N$ is the number of system states, which makes analysis of its behavior challenging. In contrast, a Lindblad-type master equation automatically ensures the physicality of the density matrix, and at the same time allows for a simpler analysis since it is characterized by a single Hamiltonian and a set of decay operators, all of dimensions $N\times N$.

§.§ Lindblad master equation with full secularization

The usual procedure to obtain a Lindblad equation from <ref> is the so-called secular approximation which consists in eliminating every term where \(\omega_{ac}\neq\omega_{bd}\). Doing so yields

\begin{align} \nonumber
\dot{\rho} = &- \frac{i}{\hbar}\left[H_\mathrm{at},\rho\right]- i \sum^{(S)}_{abd}\Big[ \Lambda_{da,db}(\omega_{bd}) \ket{a}\bra{b}, \rho \Big] \\
&+ \sum^{(S)}_{abcd}\Gamma_{ca,db}(\omega_{bd})\left(\ketbra{d}{b}\rho\ketbra{a}{c} - \frac{1}{2}\Big\{ \ket{a}\braket{c}{d}\bra{b},\rho \Big\}\right)
\end{align}

where the superscript \((S)\) in the sum indicates that only terms with \(\omega_{ac}=\omega_{bd}\) are kept. In the energy shift, this is equivalent to the condition \(\omega_{a}=\omega_{b}\) since \(\ket{c}=\ket{d}\) there. The energy shift is clearly Hermitian because it is a real and symmetric matrix. The decay term can be reexpressed by grouping the sum over transitions into sets with a given frequency \(\Omega = \omega_{ac}=\omega_{bd}\), which yields

\begin{align}
\sum_{\Omega}\sum_{\alpha\beta} \gamma_{\alpha\beta}(\Omega) \left(\sigma^{\beta}_\Omega\rho \sigma^{\alpha\dagger}_\Omega - \frac{1}{2}\Big\{ \sigma^{\alpha\dagger}_\Omega \sigma^{\beta}_\Omega,\rho \Big\}\right)
= \sum_{\Omega\epsilon} \Gamma_{\epsilon}(\Omega) \left(S^{\epsilon}_\Omega\rho S^{\epsilon\dagger}_\Omega - \frac{1}{2}\Big\{ S^{\epsilon\dagger}_\Omega S^{\epsilon}_\Omega,\rho \Big\} \right).
\end{align}

Here, Greek indices \(\alpha,\beta\) indicate spatial directions, while all dipole transitions \(d^\alpha\) with a frequency difference of \(\Omega\) are combined in the transition operators \(\sigma_\Omega^\alpha = \sum^{(\Omega)}_{ab} d^\alpha_{ab} \ket{a}\bra{b} \). The right-hand side above is obtained by diagonalizing the positive definite matrix \(\gamma_{\alpha\beta}(\Omega)=\sum_{\epsilon}M_{\alpha\epsilon}^\dagger(\Omega)\Gamma_{\epsilon}(\Omega)M_{\epsilon\beta}(\Omega)\) for each transition frequency \(\Omega\) and defining \(S^\epsilon_\Omega=\sum_\alpha M_{\epsilon\alpha}(\Omega) \sigma_\Omega^\alpha\). In this last form, it is evident that the full secularization returns a Lindblad master equation. However, the only off-diagonal terms present are the ones connecting degenerate states. This approximation has been shown to be inadequate in a variety of contexts <cit.>, since it indiscriminately removes the coupling between coherences (off-diagonal elements of \(\rho\)) and populations of non-degenerate states. Thus, in the system explored in the main text, relevant physics would be omitted within each Bohr level.

§.§ Lindblad master equation: derivation details and proof

We here show how to derive the Lindblad equation including off-diagonal terms between non-degenerate states used in the main text from the BR equation, <ref>. Instead of a full secularization as discussed above, we start by performing a partial secularization to discard terms where the timescale induced by the environment, \(\tau_E\sim\min(|\mathbf{d}^2\lambda|^{-1},|\mathbf{d}^2\gamma|^{-1})\), is much larger than that of the atomic transitions, \(\tau_\mathrm{at}\sim|\omega_{ac}-\omega_{bd}|^{-1}\). This is not an important step, but it significantly simplifies the resulting expressions. In <ref>, this is fulfilled for terms where \(\ket{a}\) and \(\ket{b}\) belong to different Bohr levels: First, if either \(\ket{c}\) or \(\ket{d}\) belongs to a different Bohr level than \(\ket{a}\) and \(\ket{b}\), respectively, then \(\tau_\mathrm{at}\) is very small compared with \(\tau_E\), and secularization is well-justified. If, instead, \(\ket{c}\) and \(\ket{d}\) belong to the same Bohr levels as \(\ket{a}\) and \(\ket{b}\), respectively, then \(\tau_E\propto 1/\gamma(\Delta_F/\hbar)\) becomes extremely large because \(\Delta_F\) is on the scale of the fine structure splitting and the spectral density approaches \(0\) when \(\omega\) goes to \(0\). Hence, even if \(\tau_\mathrm{at}\sim \hbar/\Delta_F \) is large, \(\tau_E\) is even larger in the system studied here. This partially secularized BR equation is significantly simpler than the full one, but is not yet in Lindblad form.

We next apply the geometric mean replacement discussed in the main text to <ref> and obtain

\begin{align} \nonumber
\dot{\rho} =& - \frac{i}{\hbar}\left[H_\mathrm{at},\rho\right] - i \Bigg[ \sum^{(s)}_{abcd} \tilde{\Lambda}_{ca,db} \ket{a}\braket{c}{d}\bra{b},\rho \Bigg] \\
&+ \sum^{(s)}_{abcd} \tilde{\Gamma}_{ca,db} \Bigg(\ketbra{d}{b}\rho\ketbra{c}{a} - \frac{1}{2} \Big\{  \ket{a}\braket{c}{d}\bra{b},\rho \Big\} \Bigg), \label{eq.mccauley}
\end{align}

where the superscript \((s)\) in the sum indicates the partial secularization mentioned above. This procedure is accurate for the following reason. When the spectral density, given in <ref>, is slowly varying: \(J(\omega + |\mathbf{d}^2\lambda|)\simeq J(\omega)\) and \(J(\omega + |\mathbf{d}^2\gamma|)\simeq J(\omega)\). In that case, for each term in <ref> where \(|\omega_{ac}-\omega_{bd}|< \max(|\mathbf{d}^2\gamma|,|\mathbf{d}^2\lambda|)\), the change in the value of the element is small, and the geometric mean is a good approximation. For terms where \(|\omega_{ac}-\omega_{bd}| > \max(|\mathbf{d}^2\gamma|,|\mathbf{d}^2\lambda|)\), the value might change appreciably, but its effect on the dynamics is small due to the difference in energy scales. In fact, such terms could be eliminated through an additional secularization to a good approximation.

After the replacement, the CP Hamiltonian, \(H_\textrm{CP} = \hbar \sum^{(s)}_{abc} \tilde{\Lambda}_{ca,cb} \ket{a} \bra{b} \), contains both diagonal and off-diagonal matrix elements. We note that the effect of the CR terms of the light-matter coupling is manifested in the precise values of the matrix elements, which change significantly depending on whether the CR interactions are included or not. Careful analysis shows that this Hamiltonian is Hermitian if \(\lambda(\omega_{ca})\) and \(\lambda(\omega_{bd})\) have the same sign. Since \(\lambda(\omega)\) has sign changes, in general, this is not always necessarily satisfied, leading to a potentially problematic, non-Hermitian CP Hamiltonian. In the cases studied in the manuscript, the partial secularization we performed earlier ensures that only terms with \(\omega_{ac} \simeq \omega_{bc}\) survive, and combined with the slow-varying property of the spectral density, the sign condition is satisfied always. A discussion of the general situation where this is not necessarily true will be presented in a future work.

In order for <ref> to be a Lindblad-type master equation, the decay rate tensor \(\tilde{\Gamma}_{ca,db}\) interpreted as a matrix in the combined indices \(ca\) and \(db\), sometimes called Kossakovski matrix, has to be positive semidefinite. Then, it can be diagonalized with non-negative eigenvalues and the last term in <ref> can be rewritten as a sum of standard Lindblad decay terms. While it is symmetric by construction, we are not aware of a general proof of positive semidefiniteness of the decay tensor that results when the procedure described above is applied to arbitrary spectral densities and atomic spectra. For the cases we treat in the manuscript, where the Green tensor is diagonal and cylindrically symmetric such that its Cartesian components satisfy \(G_{xx} = G_{yy}\), we give a proof below through explicit construction of the diagonalized form. Under this assumption, the decay tensor has the form \(\boldsymbol{\gamma}(\omega) = \operatorname{diag}(\gamma_{xx}(\omega), \gamma_{xx}(\omega), \gamma_{zz}(\omega))\). We now express \(\Gamma_{ca,db}(\omega)\) in terms of the spherical basis defined by

\begin{equation}
\mathbf{d}' = \begin{bmatrix} d^{+1} \\ d^{-1} \\ d^0 \end{bmatrix} = \mathbf{U}\cdot\mathbf{d} =
\begin{bmatrix}
-1/\sqrt{2} & -i/\sqrt{2} & 0 \\ 1/\sqrt{2} & -i/\sqrt{2} & 0 \\
0 & 0 & 1
\end{bmatrix}\cdot\begin{bmatrix} d^x \\ d^y \\ d^z \end{bmatrix}.
\end{equation}

By construction, the spherical components of the dipole operator, \(d^{\delta}\), connect states with a given \(m_j\) to states with \(m_j+\delta\). Due to its symmetry, \(\boldsymbol{\gamma}(\omega)\) is invariant under transformation to the spherical basis, \(\boldsymbol{\gamma}'(\omega) = \mathbf{U}\cdot\boldsymbol{\gamma}(\omega)\cdot\mathbf{U}^\dagger = \boldsymbol{\gamma}(\omega)\).
Since \(m_j\) is a well-defined quantum number of our basis states, the advantage of the spherical basis is that every transition operator \(\ketbra{a}{c}\) allowed by the selection rules (see Fig. 1b of the main text) is mediated by only one of \(d^{+1}\), \(d^{-1}\) or \(d^0\). Furthermore, because of the diagonal form of \(\boldsymbol{\gamma}'\), the transition operators \(\ketbra{d}{b}\) and \(\ketbra{c}{a}\) must have the same \(\delta\); otherwise \(\Gamma_{ca,db}(\omega) = 0\). As a consequence, we can expand the last term of <ref> as three separate sums, one for each value of \(\delta\), indicated below with the label \(\delta\) on the second summation sign:

\begin{align} \nonumber
&\sum_\delta\sum^{(s,\delta)}_{abcd} \sqrt{\Gamma_{ac,db}(\omega_{bd})}\sqrt{\Gamma_{ac,db}(\omega_{ac})}  \left[\ketbra{d}{b}\rho\ketbra{a}{c} - \frac{1}{2}\Big\{  \ket{a}\braket{c}{d}\bra{b}, \rho \Big\}\right] =\\ \nonumber
%&\sum_\delta\sum^{(s,\delta)}_{abcd} \sqrt{d_{ca}^{\delta*}\gamma_{\delta\delta}(\omega_{bd})d_{db}^\delta}\sqrt{d_{ca}^{\delta*}\gamma_{\delta\delta}(\omega_{ac})d_{db}^\delta} \left[\ketbra{d}{b}\rho\ketbra{a}{c} - \frac{1}{2}\Big\{  \ket{a}\braket{c}{d}\bra{b}, \rho \Big\} \right] = \\ \nonumber
%&\sum_\delta\sum^{(s,\delta)}_{abcd} \sqrt{d_{ca}^{\delta*}d_{db}^\delta}\sqrt{d_{ca}^{\delta*}d_{db}^\delta} \sqrt{\gamma_{\delta\delta}(\omega_{bd})}\sqrt{\gamma_{\delta\delta}(\omega_{ac})}\left[\ketbra{d}{b}\rho\ketbra{a}{c} - \frac{1}{2}\Big\{  \ket{a}\braket{c}{d}\bra{b}, \rho \Big\} \right] = \\ \nonumber
&\sum_\delta\sum^{(s,\delta)}_{abcd} d_{ca}^{\delta*}d_{db}^\delta \sqrt{\gamma_{\delta\delta}(\omega_{bd})}\sqrt{\gamma_{\delta\delta}(\omega_{ac})}\left[\ketbra{d}{b}\rho\ketbra{a}{c} - \frac{1}{2}\Big\{ \ket{a}\braket{c}{d}\bra{b}, \rho \Big\} \right] =\\
&\sum_{\delta n}\left(\Sigma^{(n)}_\delta\rho\Sigma_\delta^{(n)\dagger} -\frac{1}{2} \Big\{\Sigma_\delta^{(n)\dagger}\Sigma^{(n)}_\delta,\rho \Big\} \right),
\end{align}

where we have used that \(\gamma_{\delta\delta}(\omega) > 0\) and that the \(d^\delta_{ca}\) are real for any pair \(ca\), and in the last step, we have defined the summed transition operator

\begin{equation}
\Sigma_{\delta}^{(n)} = \sum_{db}^{(n)}  d_{db}^\delta  \sqrt{\gamma_{\delta\delta}(\omega_{bd})} \ketbra{d}{b}.
\end{equation}

Here, the states \(\ket{b}\) belong to the same Bohr level with main quantum number equal to \(n\), while the \(\ket{d}\) states can belong to any Bohr level. This simplification is a consequence of the partial secularization explained at the beginning of this subsection. In this form, the decay term is given by an explicit Lindblad operator in terms of just three decay operators for each Bohr level \(n\). We note that identical manipulations can be done on the energy shift terms, which can be refactored as

\begin{equation}
D_\delta^{(n)} = \sum_{db}^{(n)}  d_{db}^\delta  \sqrt{\lambda_{\delta\delta}(\omega_{bd})} \ketbra{d}{b}.
\end{equation}

Finally, we can rewrite <ref> in our system as

\begin{equation} \label{eq.lindblad}
\dot{\rho} = - \frac{i}{\hbar}\left[H_\mathrm{at} + H_\mathrm{CP},\rho\right] + \sum_{\delta n} L_{\Sigma^{(n)}_\delta}[\rho]
\end{equation}

where \(H_\mathrm{CP} = \hbar\sum_{\delta n}D_\delta^{(n)\dagger}D_\delta^{(n)}\) and \(L_{A}[\rho] = A \rho A^\dagger - \frac12 \{A^\dagger A, \rho\}\) is a standard Lindblad decay term. This is indeed a Lindblad equation for the atom that includes the relevant off-diagonal couplings both in the CP shift and the decay term.

§ DERIVATION OF THE EFFECTIVE HAMILTONIAN

Any Lindblad equation \(\dot{\rho} = - \frac{i}{\hbar}\left[H,\rho\right] + \sum_j L_{A_j}[\rho]\) can be rewritten as \(\dot{\rho} = - \frac{i}{\hbar} \left(\Heff\rho - \rho\Heff^\dagger\right) + \sum_j A_j \rho A_j^\dagger\), with the effective non-Hermitian Hamiltonian \(H_\mathrm{eff} = H - \frac{i}{2} \sum_j A_j^\dagger A_j \), and the terms of the last sum commonly referred to as the “refilling” or “quantum jump” terms. In physical situations where the refilling terms are negligible, the dynamics are then fully characterized by the eigenstates and eigenvalues of the effective Hamiltonian <cit.>. In the main text, we are concerned with the dynamics within a given Bohr level, in particular \(n=7\). Due to the partial secularization we performed, the effective Hamiltonian associated with <ref> is block-diagonal in Bohr levels, such that \(n\) remains a good quantum number and \([\Heff, \mathcal{P}_n] = 0 \), where \(\mathcal{P}_{n}\) is a projection operator onto the subspace with principal quantum number \(n\). Projecting the Lindblad master equation onto this subspace gives

\begin{gather}
\dot{\rho}_{n} = -\frac{i}{\hbar} \left(H^{(n)}_\mathrm{eff}\rho_{n} -\rho_n H^{(n)\dagger}_\mathrm{eff}\right) + \sum_{\delta n'}\mathcal{P}_{n}\Sigma^{(n')}_\delta\rho\Sigma_\delta^{(n')\dagger}\mathcal{P}_{n},\\
H_\mathrm{eff}^{(n)} = \mathcal{P}_n \Heff \mathcal{P}_n = H_\mathrm{at}^{(n)} + \hbar \sum_\delta \left(D_\delta^{(n)\dagger}D_\delta^{(n)}-\frac{i}{2}\Sigma_\delta^{(n)\dagger}\Sigma_\delta^{(n)}\right), \label{eq.effH}
\end{gather}

where \(\rho_n = \mathcal{P}_n \rho \mathcal{P}_n \) and \(H_\mathrm{at}^{(n)} = \mathcal{P}_n H_\mathrm{at}\mathcal{P}_n\). We thus only need to show that the refilling terms are negligible for the dynamics within a given Bohr level. To this end, we can rewrite them as

\begin{equation}
\mathcal{P}_{n}\Sigma^{(n')}_\delta\rho\Sigma_\delta^{(n')\dagger}\mathcal{P}_{n} = \sum^{(s,\delta)}_{abcd} d_{ca}^{\delta*}d_{db}^\delta \sqrt{\gamma_{\delta\delta}(\omega_{bd})}\sqrt{\gamma_{\delta\delta}(\omega_{ac})} \mathcal{P}_{n}\ket{d}\bra{b}\rho\ket{a}\bra{c} \mathcal{P}_{n}.
\end{equation}

Because of the partial secular approximation, \(\ket{a}\) and \(\ket{b}\) belong to the same Bohr level \(n'\), and due to the projection operators \(\mathcal{P}_n\), \(\ket{c}\) and \(\ket{d}\) also have the same principal quantum number, \(n\). Also, because \(\boldsymbol{\gamma}\) is only non-zero at positive frequencies, we have that \(n'\geq n\). We can immediately discard terms with \(n'>n\): They refer to the population that flows into the Bohr level \(n\) through spontaneous emission from higher-lying Bohr levels, but since we assume that the initial atomic state is in level \(n\) and there are no processes leading to higher levels, these terms do not contribute. For the remaining terms with \(n'=n\), the atomic time scales are \(\tau_\mathrm{at} \sim \hbar/\Delta_F \), but the decay-induced time scales are \(\tau_E\propto 1/\gamma(\Delta_F/\hbar)\). Given the spectral density used in the main text, in our system \(\tau_E\gg\tau_\mathrm{at}\). Thus, the effect of the terms with \(n'=n\) is negligible, and we can safely remove the “refilling” term and write the dynamics in the subspace with principal quantum number \(n\) as

\begin{equation} \label{eq.effH_evol}
\dot{\rho}_{n} = -\frac{i}{\hbar} \left(H^{(n)}_\mathrm{eff}\rho_{n} -\rho_n H^{(n)\dagger}_\mathrm{eff}\right),
\end{equation}

which is equivalent to the Schrödinger equation \(\partial_t \ket{\psi(t)} = -\frac{i}{\hbar} H_\mathrm{eff}^{(n)} \ket{\psi(t)}\).

§.§ Angular momentum conservation

For each Bohr level \(n\), the effective Hamiltonian <ref> derived above is a block diagonal matrix, with each block corresponding to a given value of the \(z\)-projection \(m_j\) of the atomic angular momentum. This is easy to see since \(H_\mathrm{at}\) conserves angular momentum, while the operators \(D_\delta^{(n)}\) and \(\Sigma_\delta^{(n)}\) connect \(m_j\) to \(m_j + \delta\), and their Hermitian conjugates connect \(m_j + \delta\) back to \(m_j\), such that overall, \(m_j\) is conserved. In contrast, physically and in the full Lindblad master equation <ref>, it is only the \(z\)-projection of the total angular momentum of the photons and atom together that is conserved due to the cylindrical symmetry of the system. Indeed, the complete master equation in <ref> does connect different \(m_j\) subspaces through the refilling term. Since we have shown this term to be negligible for the dynamics within a given subspace, we can exploit conservation of \(m_j\) to analyze its subspaces separately, and have done so in the main text by fixing \(m_j=1/2\).

§ NUMERICAL CHECK

(a) Model spectral density, \(J(\omega)\). (b) Integral of the spectral density that appears in the shift, \(\lambda(\omega)\).

In order to verify the validity of the derived Lindblad equation and effective Hamiltonian, we here apply it to a simplified system for which an exact solution is possible. To do so, we study the populations of the states with \(n=3\) in the hydrogen atom coupled to an electromagnetic bath whose spectral density is a Lorentzian. The density and the corresponding energy shift integral are shown in <ref> and <ref>, and are given by

\begin{align} \label{eq.spectral}
J(\omega) &= \frac{g^2}{\pi} \frac{\kappa/2}{(\omega-\omega_M)^2+(\kappa/2)^2}, \\
\lambda(\omega) &= \int_{-\infty}^\infty \mathrm{d}\omega' \frac{J(\omega')}{\omega-\omega'} = g^2 \frac{\omega-\omega_M}{(\omega-\omega_M)^2+(\kappa/2)^2},
\end{align}

with parameter values \(\hbar g = (9/\sqrt{5})\cdot 10^{-4}~\mathrm{eV}/(ea_0)\), \(\hbar\kappa=2\cdot 10^{-3}\) eV and \(\hbar\omega_M = 1.95\) eV. It is well-known that a Lorentzian spectral density is completely equivalent to a single mode coupled to a completely flat, i.e. Markovian, bath <cit.>, with dynamics described exactly by a Lindblad equation <cit.>,

\begin{equation} \label{eq.exact}
\dot{\rho} = -\frac{i}{\hbar}\left[H_\mathrm{at}+\hbar\omega_M a^\dagger a,\rho \right] - \frac{\kappa}{2}\{a^\dagger a,\rho\} + \kappa a\rho a^\dagger,
\end{equation}

where \(a\) is the bosonic annihilation operator of the bath mode. Hence, we can compare the approximate solutions obtained with our approaches, <ref> and <ref>, to the exact dynamics given by <ref>. We take \(\ket{\psi(0)}=\ket{n=3,l=0,j=1/2,m_j=1/2}\ket{n_\mathrm{ph}=0}\) as the initial state and propagate it in time.

Time evolution of the atomic populations. Thick brown lines: numerical solution to the exact dynamics (<ref>). Green dashed lines: effective Hamiltonian (<ref>). Black dotted lines: Lindblad master equation (<ref>).

In the exact calculations, we include the first 4 Bohr levels of the hydrogen atom with their complete fine structure (60 states), which gives converged results. The population of the \(\ket{3,l,j,1/2}\) states are plotted in <ref>. The three largest populations are perfectly well described by both our Lindblad equation and the effective, non-Hermitian Hamiltonian. It should be noted that the dynamics calculated using the BR approach are essentially the same as those obtained from the Lindblad equation and therefore not shown separately. There are two additional lines that are present only in the exact dynamics and the Lindblad equation, with populations of the order of \(10^{-3}\). These are states that become populated through the refilling terms within the \(n=3\) subspace discussed above. These are unrealistically large here because the spectral density chosen here to enable comparison with an exact result does not obey the physical constraint \(J(\omega)=0\) for \(\omega\leq 0\). In contrast, the spectral density used in the main text obeys these physical constraints and the refilling term can indeed be discarded with much less impact.

Dynamics calculated with the effective Hamiltonian. Green lines: including the contribution of the CR terms. Brown dashed lines: the rotating wave approximation has been performed on the light-matter Hamiltonian.

Having checked the validity of <ref>, we may use it to reinforce the claim that the CR terms are important in our work. In <ref>, we compare the dynamics when the effect of the CR terms is included and when it is not due to the rotating wave approximation. Clearly, the marked differences in the oscillations indicates that the CR terms significantly contribute to the CP shift and thus to the dynamics.


# Three-loop topology analysis of neutral B-meson mixing with tapir

Marvin Gerlach  Institut für Theoretische Teilchenphysik, Karlsruhe Institute
of Technology (KIT), Wolfgang-Gaede-Straße 1, 76128 Karlsruhe, Germany
<EMAIL_ADDRESS>
###### Abstract

Modern advances in particle physics depend strongly on the usage of reliable
computer programs. In this context two issues become important: The usage of
powerful algorithms to handle the amount of evaluated data properly, and a
software architecture capable to overcome the problems of maintainability and
extendability. We present our approach to such a computer program, called
tapir. This tool assists computations in perturbative quantum field theory in
many ways. Such calculations often involve the evaluation of a large amount of
Feynman diagrams with multiple loops. tapir helps in reducing the number of
diagrams, and the resulting integrals thereof, by identifying and minimizing
their topological structure. We will focus on a three-loop calculation which
is needed for the next-to-next-to leading order predictions of neutral
$B$-meson systems. We show how tapir can be utilized for this kind of
calculation.

##  1 Introduction

The outstanding experimental effort of recent years led to an unreached
precision in flavor observables. Among the most prominent ones is the mixing
of neutral meson states such as the $B_{s}$ and $B_{d}$. We take in this work
only $B_{s}$ into account, the extension to $B_{d}$ is straightforward. A more
elaborate overview for this topic is given in e.g. Ref. [1].

$B_{s}$ can oscillate in its antiparticle state $\overline{B}_{s}$ via box
diagrams. Therefore, the mass eigenstates differ from the flavor eigenstates.
A few observables for this flavor oscillation process are of particular
interest in connection to precision tests of the Standard Model. They can be
derived from the transition amplitude from one flavor eigenstate to another,
which we denote as

$\displaystyle-i(2\pi)^{4}\delta^{(4)}\left(p_{\text{i}}-p_{\text{j}}\right)\Sigma_{\text{ij}}^{s}~{}=~{}\frac{\bra{\text{i}}i\mathcal{T}\ket{\text{j}}}{2M_{B_{s}}}\,.$
(1)

i and j describe either a particle state $B_{s}$ or an antiparticle state
$\overline{B}_{s}$. $M_{B_{s}}$ is the average mass of the meson.

Especially the off-diagonal elements of this amplitude are interesting for the
mixing due to the connection to the mass and decay width matrix:

$\displaystyle M_{12}^{s}{}$
$\displaystyle=~{}\frac{\Sigma_{12}^{s}+(\Sigma_{21}^{s})^{*}}{2}~{}\equiv~{}\text{Disp}\\!\left(\Sigma_{12}^{s}\right)\,,$
(2) $\displaystyle\Gamma_{12}^{s}{}$
$\displaystyle=~{}i\left(\Sigma_{12}^{s}-(\Sigma_{21}^{s})^{*}\right)~{}\equiv~{}2\text{Abs}\\!\left(\Sigma_{12}^{s}\right)\,.$

In order to compute these off-diagonal matrix elements, it is hence necessary
to evaluate the dispersive and the absorptive part of the process
$\overline{B}_{s}\rightarrow B_{s}$. They correspond to the real and imaginary
part of the scattering amplitude where the CKM matrix elements are factored
out. According to the optical theorem, only virtual particles contribute to
$\Gamma_{12}^{s}$ which can be produced on-shell. For $M_{12}^{s}$ exists no
such restriction. Hence, $\Gamma_{12}^{s}$ is straightforward to compute when
connected to the calculation of the discontinuity of the complex plane, and
with it to Cutkosky’s cutting rules [2].

From $M_{12}^{s}$ and $\Gamma_{12}^{s}$ we are able to compute three physical
observables: The mass difference between heavy (H) and light (L) mass
eigenstates,

$\displaystyle\Delta
M_{s}~{}=~{}M_{\text{H}}^{s}-M_{\text{L}}^{s}~{}=~{}2|M_{12}^{s}|\,,$ (3)

the decay width difference,

$\displaystyle\Delta\Gamma_{s}~{}=~{}\Gamma_{\text{L}}^{s}-\Gamma_{\text{H}}^{s}~{}\approx~{}2|\Gamma_{12}^{s}|\,,$
(4)

and the CP asymmetry

$\displaystyle
a_{\text{fs}}^{s}~{}=~{}\text{Im}\left(\frac{\Gamma_{12}^{s}}{M_{12}^{s}}\right)\,.$
(5)

$a_{\text{fs}}^{s}$ is a measure of CP breaking in the mixing by quantifying
how much the mass eigenstates differ from the CP eigenstates.

$\Delta\Gamma_{s}$ is computed in an effective field theory with the effective
interaction Hamiltonian $\mathcal{H}^{\Delta B=1}$:

$\displaystyle\Gamma_{12}^{s}~{}=~{}\frac{1}{2M_{B_{s}}}\text{Abs}\\!\left(i\bra{B_{s}}\int\\!\text{d}^{4}x\
\mathcal{H}^{\Delta B=1}(x)\mathcal{H}^{\Delta
B=1}(0)\ket{\overline{B}_{s}}\right)\,.$ (6)

For $\mathcal{H}^{\Delta B=1}$ we choose the basis of Ref. [3] which is closed
under renormalization when including QCD effects. It also has the advantage of
preventing traces including $\gamma_{5}$. The numerically most relevant
effective operators are the so-called current-current operators. Taking only
those into account, we have

$\displaystyle\mathcal{H}^{\Delta
B=1}~{}=~{}\frac{4G_{F}}{\sqrt{2}}\sum_{q_{1},q_{2}}V_{q_{1}b}V_{q_{2}s}^{*}\left(C^{q_{1}q_{2}}_{1}P^{q_{1}q_{2}}_{1}+C^{q_{1}q_{2}}_{2}P^{q_{1}q_{2}}_{2}\right)\,.$
(7)

The sum iterates over $q_{1},q_{2}\in\\{u,c\\}$. The operators are given by

$\displaystyle P^{q_{1}q_{2}}_{1}{}$ $\displaystyle=~{}\left(\overline{s}\
T^{a}\gamma_{\mu}P_{L}\ q_{2}\right)\left(\overline{q_{1}}\
T^{a}\gamma^{\mu}P_{L}\ b\right)\,,$ (8) $\displaystyle P^{q_{1}q_{2}}_{2}{}$
$\displaystyle=~{}\left(\overline{s}\ \gamma_{\mu}P_{L}\
q_{2}\right)\left(\overline{q_{1}}\ \gamma^{\mu}P_{L}\ b\right)\,,$

with $P_{L}=(1-\gamma_{5})/2$.

Equation 6 can be even further simplified by using the heavy quark expansion
(HQE). This is an operator product expansion around the small parameter
$\Lambda/m_{b}$. It translates the non-local hadronic matrix element of eq. 6
to local operator insertions with new Wilson coefficients $H$ and
$\tilde{H}_{S}$. We hence get

$\displaystyle\Gamma_{12}^{s}~{}=~{}\frac{G_{F}^{2}m_{b}^{2}}{24\pi
M_{B_{s}}}\left(H\bra{B_{s}}Q\ket{\overline{B}_{s}}+\tilde{H}_{S}\bra{B_{s}}\tilde{Q}_{S}\ket{\overline{B}_{s}}\right)\,,$
(9)

with

$\displaystyle Q~{}=~{}4\left(\overline{s}_{\text{i}}\ \gamma_{\mu}P_{L}\
b_{\text{i}}\right)\left(\overline{s}_{\text{j}}\ \gamma^{\mu}P_{L}\
b_{\text{j}}\right)\,,\quad\tilde{Q}_{S}~{}=~{}4\left(\overline{s}_{\text{i}}\
P_{L}\ b_{\text{j}}\right)\left(\overline{s}_{\text{j}}\ P_{L}\
b_{\text{i}}\right)\,.$ (10)

Corrections to eq. 6 including only current-current operators are currently
known up to $\mathcal{O}(\alpha_{s})$ [4, 5, 6, 7] and only partially known to
$\mathcal{O}(\alpha_{s}^{2})$ [8, 9]. To complete the missing NNLO
contributions, diagrams as in fig. 1 have to be considered.

Figure 1: Sample diagrams which contribute at NNLO to $\Delta\Gamma_{s}$
within the effective $\Delta B=1$ theory.

##  2 Diagram and topology analysis

Analyzing the diagrams of fig. 1 gives a hint which problems may occur. In
order to allow for one-particle reducible diagrams (diagram 1), we have to use
the options “notadpoles” and “onshell” of qgraf [10] during diagram
generation. Unfortunately, “onshell” also discards diagram 2. Therefore, we
define a second diagram class where only the option “offshell” is specified.
We then use a special diagram filter of tapir [11] to remove the unwanted
diagrams. We do this by specifying the following in the  config file:

⬇

* tapir.filter external_self_energy_bridge : false

This filter excludes all diagrams whose external legs have self-energy
corrections. Flavor mixing of the self-energy, as in diagram 2, is explicitly
not effected.

qgraf generates $\mathcal{O}(4500)$ three-loop diagrams for the first diagram
class. Unfortunately, qgraf can only distinguish between interactions with
different particles, i.e. $P_{1}^{cc}$ and $P_{2}^{cc}$ cannot be
differentiated. The separation is first applied if symbolic graph
representations ( dia/ edia) are generated. Then, tapir iterates over all
possible interaction combinations as stated in the  vrtx file. Thus, the
number of diagrams increases to $\mathcal{O}(17000)$ if all current-current
operators of eq. 8 are regarded. The second class starts with
$\mathcal{O}(3000)$ diagrams generated by qgraf, but reduces to
$\mathcal{O}(2000)$ when applying the described diagram filter. The actual
number of remaining diagrams of the second class with distinguished operators
is $\mathcal{O}(8000)$.

A further simplification concerns the specific kinematics of the process
$\overline{b}(q_{1})s(q_{2})\rightarrow\overline{s}(q_{3})b(q_{4})$. Since we
match the results of two effective field theories to get $H$ and
$\tilde{H}_{S}$, the result must not depend on the kinematics of the external
particles. Thus, we can choose the three-momenta of the $s$-quarks as
$\vec{q}_{2}=\vec{q}_{3}=0$. In the limit $m_{s}\rightarrow 0$, this kinematic
is topologically equal to diagrams where the external $s$-lines are removed.
We can make use of this in tapir with

⬇

* tapir.external_momentum q2:0

* tapir.external_momentum q3:0

The external momenta are labeled according to the order of how the external
particles were declared in  qgraf.dat. Additionally, we specify that the $b$\-
and $c$-quarks are massive:

⬇

* tapir.mass fb:M1

* tapir.mass fc:M2

After specifying a  prop and  vrtx file with all symbolic Feynman rules, we
analyze and minimize the occurring topologies. We call tapir with the
following command line options:

⬇

$ /path/to/tapir -c class1.conf -q qlist.3 -m -t class1.top -k 8

The “-c” argument specifies the  config file. The argument of “-q” defines the
qgraf output file that was generated using the style file which comes with
tapir. “-m” minimizes the topologies of the diagrams and tries to map them
onto each other. The option “-t” outputs the remaining topologies as a  topsel
file. This file format was adopted from q2e and exp [12, 13, 14] and serves as
the standard format for graph topologies around tapir. The last command line
option “-k” demands tapir to use 8 cores for parallel evaluation. The
execution per diagram class takes $\mathcal{O}($10\text{\,}\mathrm{s}$)$ on a
modern desktop PC.

##  3 Topology minimization

The implemented minimization starts with extracting the topological features
of each diagram. This includes the connections of the underlying graph, i.e.
which vertices are connected by which edge. It also takes the edge coloring
into account, i.e. the mass of the propagating particle. This information is
then encoded in a label which is unique for a given Feynman graph topology:
the Nickel index [15]. With this index at hand it is straightforward to
identify a priori different graphs, by using it as a hash function and
comparing the generated hashes. This procedure enables a minimization in a
highly parallelized manner.

The two diagram classes defined above lead minimized to $\mathcal{O}(1000)$
and $\mathcal{O}(400)$ remaining topologies, respectively. In the next step we
want to combine both topology sets, and reduce them even further. For this
purpose, we create a new  config file with the following entries:

⬇

* tapir.topo_ignore_bridges

* tapir.topo_remove_duplicate_lines

Here, we use a set of destructive topology filters, i.e. they change the graph
topology, but do not affect the underlying integral families. The option
“tapir.topo_ignore_bridges” removes lines from a topology which do not carry a
loop momentum, so-called bridges. “tapir.topo_remove_duplicate_lines” removes
lines which have the exact same momentum and mass as another line.

As a next step, we want to use the topological information to build symbolic
integral families. For this, we generate a FORM [16] file that combines the
scalar loop-dependent factors to a topology function. This  topology file can
be embedded in a FORM-based setup, as it is common for multi-loop
calculations. Before including this file, some kind of tensor reduction has to
be applied to the symbolic expressions of the diagrams generated by tapir. The
generation of the  topology file can be varied by several options. For our
case, we may use

⬇

* tapir.topo_complete_momentum_products

to replace numerator functions of the integral family by a standard propagator
function. For example, it may happen that the product $q.k$ cannot be
expressed through the remaining propagators. Thus, $q.k$ would be part of the
integral family. The presented option replaces the product with

$q.k~{}=~{}\frac{1}{2}\left((q+k)^{2}-q^{2}-k^{2}\right)\,.$ (11)

The combination $q+k$ is then replaced by an ancillary momentum $P_{1}$, and
$1/(-P_{1}^{2})^{n}$ is included in the integral family. The minus sign is
used for convenience.

We finally can call tapir once again with

⬇

$ /path/to/tapir -c combine.conf -i class12.top -m -t final.top -f include/ -k
8 -pf

This time, we read a concatenation of the previously created  topsel files.
One has to make sure that the topologies of the different classes can be
distinguished using the option “tapir.topology_name” beforehand. The option
“-f” creates a folder with the  topology files in addition to a  topology list
file containing all integral family definitions in a Mathematica readable
format. The latter can be used as input for reduction programs like FIRE [17].

Additionally, we use the option “-pf” which usually performs a partial
fraction decomposition of the final integral families. Thus, we ensure to get
only integral families with linear independent denominator functions. This is
a necessary pre-condition for integration-by-parts algorithms. The combination
of both diagram classes with the mentioned topology filters results in
$\mathcal{O}(900)$ topologies. It turns out, the partial fraction
decomposition does not apply since the options topo_ignore_bridges and
topo_remove_duplicate_lines already remove all occurring linear dependencies
in the integral families. Nevertheless, using “-pf” is always a good
precaution.

The remaining $\mathcal{O}(900)$ topologies are the result of only a coarse
minimization. A program like exp can map the diagrams on even less topologies.
Thus, our naive minimization procedure can be regarded as a relatively fast
“pre-filtering”.

The generated  topology files, e.g. for the topology of diagram 4 in fig. 1,
have the following form:

⬇

* Reducible numerator momentum replacements

id p3 = -p4 - q1;

id p1 = p4 - p5 + p6;

id p2 = p4 - p5 + p6 + q1;

id P1 = p5 + p6;

id P2 = p4 + p6;

id P3 = p5 + q1;

.sort

* Numerator momentum product replacements

id p6.q1 = -p1.p1/2 + p2.p2/2 - p3.p3/2 + p4.p4/2 + p5.q1;

id p4.q1 = p3.p3/2 - p4.p4/2 - q1.q1/2;

id p4.p5 = -p1.p1/2 + p4.p4/2 + p4.p6 + p5.p5/2 - p5.p6 + p6.p6/2;

id p5.p6 = (P1.P1 - p5.p5 - p6.p6)/2;

id p4.p6 = (P2.P2 - p4.p4 - p6.p6)/2;

id p5.q1 = (P3.P3 - p5.p5 - q1.q1)/2;

.sort

* Define massive propagators

id p2.p2 = -1/s2m1 + M1^2;

id p4.p4 = -1/s4m1 + M1^2;

id p5.p5 = -1/s5m2 + M2^2;

id p6.p6 = -1/s6m2 + M2^2;

.sort

* Combine to scalar topology function

id s6m2^n0? * s5m2^n1? * s4m1^n2? * s2m1^n3? * 1/p1.p1^n4? * 1/p3.p3^n5? *
1/P1.P1^n6? * 1/P2.P2^n7? * 1/P3.P3^n8? =

(-1)^n4 * (-1)^n5 * (-1)^n6 * (-1)^n7 * (-1)^n8 *

DB1TopoNNLOdr3l4380(n0,n1,n2,n3,n4,n5,n6,n7,n8);

.sort

Herewith a tensor reduced, symbolic diagram is straightforwardly replaced by
an integral family in FORM. The corresponding  topology list entry is given by

⬇

{"DB1TopoNNLOdr3l4380", {M2^2 - p6^2, M2^2 - p5^2, M1^2 - p4^2, M1^2 - (p4 -
p5 + p6 + q1)^2, -(p4 - p5 + p6)^2, -(p4 + q1)^2, -(p5 + p6)^2, -(p4 + p6)^2,
-(p5 + q1)^2}, {p4, p5, p6}}

which is in the form of

⬇

{"Family name", List of denominators, List of loop momenta}

##  4 Conclusion

We have shown how the program tapir can be used to simplify Feynman graph
topologies and how they can be expressed as symbolic integral families. We
have considered the example of $B_{s}-\overline{B}_{s}$ mixing at three-loop
order. The program was already used in similar two-loop calculations [18, 19].

The program and an extensive documentation are freely available at the gitlab
page: https://gitlab.com/F.Herren/tapir.

I thank Matthias Steinhauser for careful proofreading of this article. I also
want to thank Florian Herren and Martin Lang for the good collaborative work
on tapir.

This research was supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under grant 396021762 — TRR 257 “Particle Physics
Phenomenology after the Higgs Discovery”.

## References

## References

  * [1] Buras A J 2020 Gauge Theory of Weak Decays: The Standard Model and the Expedition to New Physics Summits (Cambridge University Press)
  * [2] Cutkosky R E 1960 Journal of Mathematical Physics 1 429–433
  * [3] Chetyrkin K G, Misiak M and Munz M 1998 Nucl. Phys. B 520 279–297 (Preprint hep-ph/9711280)
  * [4] Beneke M, Buchalla G, Greub C, Lenz A and Nierste U 1999 Phys. Lett. B 459 631–640 (Preprint hep-ph/9808385)
  * [5] Ciuchini M, Franco E, Lubicz V, Mescia F and Tarantino C 2003 JHEP 08 031 (Preprint hep-ph/0308029)
  * [6] Beneke M, Buchalla G, Lenz A and Nierste U 2003 Phys. Lett. B 576 173–183 (Preprint hep-ph/0307344)
  * [7] Lenz A and Nierste U 2007 JHEP 06 072 (Preprint hep-ph/0612167)
  * [8] Asatrian H M, Hovhannisyan A, Nierste U and Yeghiazaryan A 2017 JHEP 10 191 (Preprint hep-ph/1709.02160)
  * [9] Asatrian H M, Asatryan H H, Hovhannisyan A, Nierste U, Tumasyan S and Yeghiazaryan A 2020 Phys. Rev. D 102 033007 (Preprint hep-ph/2006.13227)
  * [10] Nogueira P 1993 J. Comput. Phys. 105 279–289
  * [11] Gerlach M, Herren F and Lang M 2022 (Preprint hep-ph/2201.05618)
  * [12] Harlander R, Seidensticker T and Steinhauser M 1998 Phys. Lett. B 426 125–132 (Preprint hep-ph/9712228)
  * [13] Seidensticker T 1999 Automatic application of successive asymptotic expansions of Feynman diagrams 6th International Workshop on New Computing Techniques in Physics Research: Software Engineering, Artificial Intelligence Neural Nets, Genetic Algorithms, Symbolic Algebra, Automatic Calculation (Preprint hep-ph/9905298)
  * [14] Harlander R, Seidensticker T and Steinhauser M 1997 q2e/exp URL http://sfb-tr9.ttp.kit.edu/software/html/q2eexp.html
  * [15] Nickel B, Meiron D and Baker G A J 1977 University of Guelph Report
  * [16] Ruijl B, Ueda T and Vermaseren J 2017 (Preprint hep-ph/1707.06453)
  * [17] Smirnov A V and Chukharev F S 2020 Comput. Phys. Commun. 247 106877 (Preprint hep-ph/1901.07808)
  * [18] Gerlach M, Nierste U, Shtabovenko V and Steinhauser M 2021 JHEP 07 043 (Preprint hep-ph/2106.05979)
  * [19] Gerlach M, Nierste U, Shtabovenko V and Steinhauser M 2022 JHEP 04 006 (Preprint hep-ph/2202.12305)