# Nonlocality of the energy density of a spontaneously emitted single-photon from a Hydrogen atom

Maxime Federico, H R Jauslin

► To cite this version:

Maxime Federico, H R Jauslin. Nonlocality of the energy density of a spontaneously emitted single-photon from a Hydrogen atom. Journal of Physics A: Mathematical and Theoretical, 2024, 58 (1), pp.015304. 10.1088/1751-8121/ad97fa. hal-04906021

HAL Id: hal-04906021

https://hal.science/hal-04906021v1

Submitted on 22 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Nonlocality of the energy density of a spontaneously emitted single-photon from a
Hydrogen atom

M. Federico¹,² and H.R. Jauslin²,∗

¹Telecom Paris, Institut Polytechnique de Paris, 19 Place Marguerite Perey, 91120 Palaiseau, France

²Laboratoire Interdisciplinaire Carnot de Bourgogne - ICB,
Université de Bourgogne, CNRS UMR 6303, BP 47870, 21078 Dijon, France

(Dated: January 9, 2025)

We analyze through the expectation value of the energy density the spatial nonlocality of single photons emitted by the spontaneous decay of a Hydrogen atom. By using a minimal coupling between the quantized electromagnetic field and the atom, we compute the state of the photon under the assumption that only a single-photon is produced. The calculations are thus performed in the subspace of single-photon states which is essentially equivalent to the rotating wave approximation. We obtain a characterization of the spatial decay of the energy density. We compute the asymptotic limit of large distances from the atom at each given time, and find an algebraic behavior of $1/r^6$. This result confirms that the energy density of single-photon states is nonlocal and the algebraic decay is far from the maximal quasiexponential localization predicted by the theory.

# I. INTRODUCTION

The description of the spontaneous decay of a Hydrogen atom is a standard calculation that can be done using either the Weisskopf-Wigner theory [1-3] or time-dependent perturbation theory, e.g., with the Fermi golden rule [3]. However, these approaches in general consider a dipolar coupling between the field and the atom, taken as a two-level system, which prevents from computing the state of the emitted photons since there is a frequency Lamb shift given by a diverging intergral. The treatment of this divergence could be done by a renormalization procedure, as suggested in [4]. This problem was discussed, e.g., in [5-10] and a way out was found by considering a minimal coupling instead of the dipole approximation, which provides an intrinsic ultraviolet cutoff function that avoids the need of a renormalization. In this article we analyze the spatial properties of the single-photon states that would be produced in this process. In view of the intrinsic spatial nonlocality of single-photons, as discussed, e.g., in [11-18], we would like to characterize this nonlocality from the perspective of spontaneous emission. To do so, we will consider the approach that uses the energy density to probe the nonlocality of photons, as it was originally done in [16, 17] and recently extended to all single-photon states in [18]. In this article it was shown that for any single-photon state, the expectation value of the energy density in a small volume is different from zero anywhere in space. This result leads to several natural questions if we consider the photons produced by some concrete source: (i) If one assumes that what is produced are single-photon states, what would be their spatial localization properties? (ii) How does the mean value of the energy density depend on the distance between the source and the detector? We analyze these questions with the example of the spontaneous emission of a Hydrogen atom prepared in its first excited state. We start with the complete Hamiltonian of the atom interacting with the quantized electromagnetic field. The hypothesis that only a single-photon is produced is implemented by projecting the model to the single-photon subspace. Furthermore we keep only the fundamental and the first excited state of the atom, corresponding to the Lyman-α transition. The time evolution in the Schrödinger picture can be represented in terms of a system of Weisskopf-Wigner equations that can be solved exactly in terms of some integrals [2]. This provides the spatio-temporal description of the photon state, which we analyze at each fixed time as a function of the distance from the source. The expectation value of the electromagnetic energy density is non-zero in any finite volume at arbitrary distances *r* from the source, as it was proven for all single-photon states in [18]. The asymptotic analysis for large distances leads to the result that at each fixed time the energy density is proportional to $1/r^6$.

We conclude that the spontaneous emission cannot produce a strictly single-photon state since it would contradict causality [19, 20]: At any fixed time *t* there would be a nonzero probability of detecting the photon with an instrument located at an arbitrary distance *r* from the source, in particular one could take *r* > *ct*, where *c* is the speed of light [21-23]. We remark that we made the analysis for the situation where the atom is in the first excited state at *t* = 0, but a similar argument can be made if at *t* = 0 the atom is in a linear combination of the ground state and the first excited state. This could correspond to the state prepared with a π-pulse of finite duration.

<footnote>∗ jauslin@u-bourgogne.fr</footnote>
Several works have already addressed related questions in the context of spontaneous emission [13, 24–28]. In particular, causality properties have been analyzed in [29–32], where different regimes, in the near- mid- and far-field regions were described in terms of wave functions constructed by Glauber’s extraction principle from a first order correlation function. This wave function is identical, up to a global constant, with the Bialynicki-Birula representation of single-photon states [33], which is in turn equivalent to the representation of the states that we will use in this work (Landau-Peierls) [34]. However, we note that the formulas (25) in Ref. [31] do not provide the asymptotic behavior for large distances $r = |\vec{x}|$, since they still contain integrals over $k$ involving functions of $|\vec{x}|$ that will modify the asymptotic form for large $|\vec{x}|$. Moreover, we stress that our spatial analysis will not be performed for the state of the field but for its energy density expectation value. This is a major difference since it has been shown (see, e.g., [18]) that in some particular situations, the state function can be spatially localized even though the corresponding energy density features a nonlocality. The nonlocality property is thus not a property of the states alone but of the joint representation of the states and the corresponding local observables.

Causality properties were also analyzed in [35] where it was shown using the electric energy density observable that the nonresonant terms are needed for the spontaneous emission to be causal. This result is in agreement with the characterization we perform in this work for the resonant part, for which we also provide an expilcit asymptotic behavior. Such an asymptotic analysis, performed through perturbation theory calculations, was also done in [36] but for the dressed vacuum only, i.e., atom in its ground state and zero photons.

Finally, as it was noted by Shirokov [37], since the eigenfunctions of the Hydrogen atom decrease exponentially in the radial variable, they are nonzero everywhere. Therefore the emision of photons cannot be expected to be local, but should have at least an exponential tail. This effect has been described as a blurring of the causal light cone. However, this is an exponential decay which is a relatively small effect compared to the nonlocality properties described below.

To analyze the spatial properties of the photon states, we will use throughout this work a position space formulation of the quantum theory of the electromagnetic field called the Landau-Peierls representation (see Appendix A for its definition). Any other position space representation could in principle be used and the result is independent of this choice [34].

The structure and the main results of the article are as follows. In Section II we formulate the Weisskopf-Wigner model with the $\vec{A} \cdot \vec{p}$ interaction and we use the vector spherical harmonics and helicity basis for single-photon states, as proposed in [5–10]. This basis has the advantage that it leads easily to exact selection rules for the interaction with Hydrogen atoms. In Section III we write the Weisskopf-Wigner solution in terms of integrals. This solution is exact under the condition that the spectrum of the associated Friedrichs operator is absolutely continuous, which is expected in the small coupling regime. Some of these integrals can be calculated explicitly if one makes a resonant approximation which has essentially the same effect as the Weisskopf-Wigner approximation, but without extending an integral to non-physical negative frequencies. With the coupled atom-electromagnetic field state, we calculate in Section IV the expectation value of the energy density at an arbitrary position in $\mathbb{R}^3$. We then analyze its dependence on the radial distance, in particular for the far-field, which characterizes the global extension of the nonlocality. We obtain the result that the expectation value of the energy density is proportional to $1/r^6$ for large $r$, i.e., it has an algebraic decay, which is far from the theoretically maximal quasiexponential decay predicted in [16].

## II. WEISSKOPF-WIGNER MODEL

We consider a non-relativistic Hydrogen atom represented by a two-level system corresponding to the Lyman-α transition, i.e., with ground state $|n_g = 1; j_g = 0; m_g = 0\rangle \equiv |g\rangle$ and excited state $|n_e = 2; j_e = 1; m_e = 0, \pm 1\rangle \equiv |e\rangle$. This transition is relevant since it can be prepared experimentally by applying, e.g., a resonant π-pulse on the ground state with a linearly polarized laser light for $m_e = 0$ and a circularly polarized laser light for $m_e = \pm 1$. The wavefunctions associated to the ground and excited states are, respectively,

$$|g\rangle = \varphi_g(r, \vartheta, \varphi) = \frac{1}{\sqrt{\pi r_B^3}} e^{-\frac{r}{r_B}}, \qquad |e\rangle = \varphi_e(r, \vartheta, \varphi) = \beta_{m_e}(\vartheta, \varphi) \frac{r}{r_B} e^{-\frac{r}{2r_B}}, \qquad (1)$$

where $r_B = \hbar/(amc)$ is the Bohr radius expressed, here, in terms of the fine structure constant $\alpha$, the electron mass $m$, and the speed of light $c$. The functions $\beta_{m_e}(\vartheta, \varphi)$ depend on the choice that is made for $m_e = \{0, \pm 1\}$ and read

$$\beta_0(\vartheta, \varphi) = \frac{1}{4\sqrt{2\pi r_B^3}} \cos \vartheta, \qquad \beta_{\pm 1}(\vartheta, \varphi) = \mp \frac{1}{8\sqrt{\pi r_B^3}} \sin \vartheta e^{\pm i\varphi}. \qquad (2)$$

The free Hamiltonian of the two-level atom can be expressed as

$$\hat{H}_{\text{at}} = E_g |g\rangle \langle g| + E_e |e\rangle \langle e|, \qquad (3)$$
where $E_g$ and $E_e$ are the energies of the ground and excited states, respectively. To describe the spontaneous emission, the two-level Hydrogen atom is interacting with the quantized electromagnetic field which we write using the basis $|k, J, M, \lambda\rangle \equiv \vec{\psi}_{k, J, M}^{(\lambda)}$ of helicity vector spherical harmonics (see Appendix B for more details about these functions)

$$ \hat{H}_{\text{em}} = \int_{0}^{\infty} dk \sum_{J,M} \sum_{\lambda=\pm} \hbar \omega_k \hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}}^{\dagger} \hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}}. \quad (4) $$

We remark, here, that the creation-annihilation operators $(\hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}}^{\dagger}, \hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}})$ are defined on functions in the Landau- Peierls representation [34, 38] (see Appendix A), and fulfill the following commutation relation $[\hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}}, \hat{B}_{\vec{\psi}_{k',J',M'}^{(\lambda')}}] = \delta(k-k')\delta_{J,J'}\delta_{M,M'}\delta_{\lambda,\lambda'}$. For the interaction part, and following [5], we take an operator of the form

$$ \hat{H}_{\text{int}} = -\frac{e}{m} \vec{p} \cdot \vec{A}(\vec{x}), \quad (5) $$

where $\vec{p} = -i\hbar\nabla$ is the momentum operator of the electron and $\vec{A}(\vec{x})$ is the vector potential operator in the Coulomb gauge. We have neglected the $\hat{A}^2$ term coming from the minimal coupling since it has been shown by Moses [5] (see also [8, 9]) that taking an interaction of the form $\vec{p} \cdot \vec{A}$ is enough to have a consistent model allowing to compute the photon state. Moreover, the typical order of magnitude for the effects related to the $\hat{A}^2$ term have been estimated to be much smaller than those from the $\vec{p} \cdot \vec{A}$ term [39]. The resulting total Hamiltonian is finally

$$ \hat{H} = \hat{H}_{\text{at}} + \hat{H}_{\text{em}} + \hat{H}_{\text{int}} \quad (6a) $$

$$ = E_g |g\rangle \langle g| + E_e |e\rangle \langle e| + \int_0^\infty dk \sum_{J,M} \sum_{\lambda=\pm} \hbar \omega_k \hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}}^{\dagger} \hat{B}_{\vec{\psi}_{k,J,M}^{(\lambda)}} - \frac{e}{m} \vec{p} \cdot \vec{A}(\vec{x}). \quad (6b) $$

In the following, since we are interested in the characterization of single photons, we project the system into the subspace generated by the vacuum and the single-photon states, which corresponds in the tensor product Hilbert space to the subspace generated by

$$ \{|e; \emptyset\rangle, |g; \vec{\psi}_{k,J,M}^{(\lambda)}\rangle\}, $$

i.e., the atom in its excited state and zero photons and the atom in its ground state and one photon. This approximation is essentially equivalent to the rotating wave approximation that is performed in general for the dipole coupling. Within this subspace, the interaction coefficient takes the form [5]

$$ \rho(k, J, M, \lambda) = -\frac{e}{m} \langle g; \vec{\psi}_{k, J, M}^{(\lambda)} | \vec{p} \cdot \vec{A}(\vec{x}) | e; \emptyset \rangle \quad (7a) $$

$$ = \left(\frac{2}{3}\right)^{7/2} \sqrt{\frac{\alpha^5}{\pi}} mc^2 \frac{k/K}{\sqrt{k}} \frac{\delta_{J,1} \delta_{M, m_e}}{\left[1 + \left(\frac{k}{K}\right)^2\right]^2} \quad (7b) $$

$$ \equiv \rho(k) \delta_{J,1} \delta_{M, m_e}, \quad (7c) $$

where $K = 3/(2r_B)$. Here, we have chosen a slightly different convention for the global phases of the eigenfunctions $\vec{\psi}_{k, J, M}^{(\lambda)}$ than the one of [5]. Our convention has the advantage that $\rho$ is real and independent of the helicity $\lambda$.

### III. SOLUTION IN TERMS OF INTEGRALS

Starting with a general state in the subspace of interest

$$ |\Psi\rangle = c_0 |e; \emptyset\rangle + \int_0^\infty dk \sum_{J,M,\lambda} c_{k,J,M,\lambda} |g; \vec{\psi}_{k,J,M}^{(\lambda)}\rangle, \quad (8) $$

we want to solve the Schrödinger equation

$$ i\hbar \frac{\partial}{\partial t} |\Psi\rangle = \hat{H} |\Psi\rangle, \quad (9) $$
with the initial condition $c_0(t = 0) = 1$ and $c_{k,J,M,\lambda}(t = 0) = 0$, for any $k, J, M$ and $\lambda$. Expressed for the nonzero coefficients, we arrive at the set of equations

$$i\hbar\dot{c}_0(t) = E_e c_0(t) + \int_0^\infty dk \sum_{\lambda=\pm} \rho(k) c_{k,\lambda}(t), \quad (10a)$$

$$i\hbar\dot{c}_{k,\lambda}(t) = \rho(k)c_0(t) + (E_g + \hbar\omega_k)c_{k,\lambda}(t), \quad (10b)$$

where we have introduced the shorter notation $c_{k,J=1,M=m_e,\lambda} \equiv c_{k,\lambda}$. Moreover, since $\rho(k)$ is $\lambda$-independent, the last equation is also $\lambda$-independent and thus $c_{k,+}(t) = c_{k,-}(t)$ for all $t$ and one can define $d_k(t) = \sqrt{2}c_{k,\lambda}(t)$, to rewrite the set of dynamical equations as

$$i\dot{c}_0(t) = \omega_a c_0(t) + \int_0^\infty dk \frac{\sqrt{2}}{\hbar} \rho(k) d_k(t), \quad (11a)$$

$$i\dot{d}_k(t) = \frac{\sqrt{2}}{\hbar} \rho(k) c_0(t) + \omega_k d_k(t), \quad (11b)$$

where we have set the energy of the ground state to $E_g = 0$, and introduced the Bohr frequency of the two-level system $\omega_a = (E_e - E_g)/\hbar$. Expressed in terms of $\omega_k = ck$, equations (11) can be written in a matrix form

$$i \begin{bmatrix} \dot{c}_0 \\ \dot{D}_{\omega_k} \end{bmatrix} = \begin{bmatrix} \omega_a & \langle \tilde{\rho} | \cdot | \cdot \rangle \\ \tilde{\rho} & \omega_k \end{bmatrix} \begin{bmatrix} c_0 \\ D_{\omega_k} \end{bmatrix}, \quad (12)$$

where $D_{\omega_k} = d_{\omega_k}/\sqrt{c}$ and

$$\tilde{\rho}(\omega_k) = \sqrt{\frac{2}{c\hbar^2}} \rho(\omega_k). \quad (13)$$

This formulation is known as the Friedrichs-Lee model [40, 41] and the exact solution in terms of integrals can be formulated as [2]

$$c_0(t) = \int_0^\infty d\omega \, g(\omega) e^{-i\omega t}, \quad (14a)$$

$$D_{\omega_k}(t) = -ie^{-i\omega_k t} \tilde{\rho}(\omega_k) \int_0^t dt' e^{i\omega_k t'} c_0(t'), \quad (14b)$$

with

$$g(\omega) = \frac{\Gamma(\omega)}{2\pi (\omega - \omega_e - \Delta(\omega))^2 + \frac{\Gamma(\omega)^2}{4}}, \quad (15)$$

and

$$\Gamma(\omega) = 2\pi|\tilde{\rho}(\omega)|^2, \quad \Delta(\omega) = \text{pv} \int_0^\infty d\omega_k \frac{|\tilde{\rho}(\omega_k)|^2}{\omega - \omega_k}, \quad (16)$$

where pv is the principal value. In the weak coupling regime, $g(\omega)$ has a high peak centered at $\omega = \omega_a$ and one can thus make the following approximation [2]

$$g(\omega) \approx g_w(\omega) = \frac{\Gamma(\omega_a)}{2\pi (\omega - \omega_a - \Delta(\omega_a))^2 + \frac{\Gamma(\omega_a)^2}{4}} \quad (17a)$$

$$= \frac{\Gamma_a}{2\pi (\omega - \omega_a - \Delta_a)^2 + \frac{\Gamma_a^2}{4}}, \quad (17b)$$

where $\Gamma_a = \Gamma(\omega_a)$ and $\Delta_a = \Delta(\omega_a)$. This approximation allows to express $c_0(t)$ as

$$c_0(t) \approx \int_0^\infty d\omega \, g_w(\omega) e^{-i\omega t} = e^{-\frac{\Gamma_a}{2}t} e^{-i(\omega_a + \Delta_a)t}, \quad (18)$$

yielding a result that is essentially equivalent to the Weisskopf-Wigner approximation. From this, one can further compute

$$D_{\omega_k}(t) \approx -ie^{-i\omega_k t} \tilde{\rho}(\omega_k) \int_0^t dt' e^{-\frac{\Gamma_a}{2}t'} e^{-i(\omega_a - \omega_k + \Delta_a)t'}. \quad (19)$$
IV. ENERGY DENSITY EXPECTATION VALUE

Following [16, 18], we want to characterize the localization of the emitted single-photon state using the energy density operator

$$ \hat{\mathcal{E}}_{\text{em}}(\vec{x}) = \frac{\varepsilon_0}{2} : (\hat{\vec{E}}^2(\vec{x}) + c^2 \hat{\vec{B}}^2(\vec{x})) :; \qquad (20) $$

for which the expectation value for a single-photon state $|\vec{\psi}\rangle$ in the Landau-Peierls (LP) representation (see Appendix A) reads [16, 18]

$$ \langle \hat{\mathcal{E}}_{\text{em}}(\vec{x}) \rangle = \hbar |\Omega^{1/2} \vec{\psi}^{(h+)}(\vec{x})|^2 + \hbar |\Omega^{1/2} \vec{\psi}^{(h-)}(\vec{x})|^2. \qquad (21) $$

The exponents $(h\pm)$ refer here to the positive and negative helicity parts of the LP representation (see Appendix A). This splitting into helicity components is an important feature of the result since it leads to the nonlocality of the energy density for all single-photon states [18] that we want to characterize for the spontaneous emission.

Using the results of the preceding section, we are going to express the state of the photon that is spontaneously emitted by the Lyman-α transition of the Hydrogen atom in the single-photon approximation, and then compute the associated expectation value of the electromagnetic energy density. We consider here only the photon part of the state, i.e.,

$$ \vec{\psi}_{\text{ph}}(r, \vartheta, \varphi) = \int_{0}^{\infty} dk \sum_{\lambda=\pm} c_{k,\lambda}(t) \vec{\psi}_{k}^{(\lambda)}(r, \vartheta, \varphi), \qquad (22) $$

where we use the abridged notation $\vec{\psi}_k^{(\lambda)} = \vec{\psi}_{k, J=1, M=m_e}^{(\lambda)}$. The positive and negative helicity parts of the state can be written as

$$ \vec{\psi}_{\text{ph}}^{(h\pm)}(r, \vartheta, \varphi) = \int_{0}^{\infty} dk \frac{d_k(t)}{\sqrt{2}} \vec{\psi}_{k}^{(\pm)}(r, \vartheta, \varphi). \qquad (23) $$

Therefore, to compute the mean value, one needs to express

$$ |\Omega^{1/2} \vec{\psi}_{\text{ph}}^{(h\lambda)}(r, \vartheta, \varphi)|^2 = \left| \int_{0}^{\infty} dk \sqrt{\frac{\omega_k}{2}} d_k(t) \vec{\psi}_k^{(\lambda)}(r, \vartheta, \varphi) \right|^2 \qquad (24) $$

for both helicities. To do so, we use the expression of $\vec{\psi}_k^{(\lambda)}$ in terms of vector spherical harmonics and Bessel functions, as given in Appendix B, to obtain

$$ |\Omega^{1/2} \vec{\psi}_{\text{ph}}^{(h\lambda)}|^2 = \sum_{L=0}^{2} F_L(r,t) \vec{Y}_{1,m_e}^{L}(\vartheta, \varphi), \qquad (25) $$

where the $F_L$ functions are

$$ F_0(r,t) = i \sqrt{\frac{c}{3\pi}} \int_0^\infty dk \, k^{3/2} d_k(t) j_0(kr), \qquad (26a) $$

$$ F_1(r,t) = \lambda \sqrt{\frac{c}{2\pi}} \int_0^\infty dk \, k^{3/2} d_k(t) j_1(kr), \qquad (26b) $$

$$ F_2(r,t) = -i \sqrt{\frac{c}{6\pi}} \int_0^\infty dk \, k^{3/2} d_k(t) j_2(kr). \qquad (26c) $$

In order to simplify the expression of the mean value, we perform an asymptotic analysis of the $F_L$ functions for distances $r$ far from the atom. The details of this analysis are given in the Appendix C. We obtain that the leading term is $F_1$ which behaves as

$$ F_1(r,t) \underset{r \to \infty}{\sim} -i\lambda \left(\frac{2}{3}\right)^{11/2} \sqrt{\frac{\alpha^5 mc^2 r_B^2}{c \pi \hbar}} \mathcal{T}(t) \frac{1}{r^3}, \qquad (27) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Plot of the angular distribution $\gamma_{|m_e|}(\vartheta)$ of the photon state emitted by spontaneous emission of the Lyman-$\alpha$ transition of a Hydrogen atom, obtained in the asymptotic limit of large distances. The green solid line is the distribution for an excited state with $m_e = 0$. The red dashed line is the distribution for an excited state with $|m_e| = 1$.</figcaption></figure>

where $\mathcal{T}(t)$ is a dimensionless function of time given in equation (C16). The expectation value of the energy density becomes thus in the asymptotic limit of large $r$

$$ \left| \Omega^{1/2} \vec{\psi}_{\text{ph}}^{(h\lambda)}(\vec{x}) \right|^2 \underset{r \to \infty}{\sim} |F_1(r,t)|^2 \vec{Y}_{1,m_e}^{1\star} \cdot \vec{Y}_{1,m_e}^1, \qquad (28) $$

which can be computed for different values of $m_e$. For $m_e = 0$, the angular part reads

$$ \vec{Y}_{1,0}^{1\star} \cdot \vec{Y}_{1,0}^1 = \frac{3}{8\pi} \sin^2 \vartheta, \qquad (29) $$

and for $|m_e| = 1$

$$ \vec{Y}_{1,\pm 1}^{1\star} \cdot \vec{Y}_{1,\pm 1}^1 = \frac{3}{16\pi} (1 + \cos^2 \vartheta). \qquad (30) $$

In the end, the expectation value as a function of $m_e$ can be expressed as

$$ \langle \hat{\mathcal{E}}_{\text{em}} \rangle = \left(\frac{2}{3}\right)^{10} \frac{m^2 c^3 r_{\text{B}}^4 \alpha^5}{2\pi^3 \hbar} \gamma_{|m_e|}(\vartheta) |\mathcal{T}(t)|^2 \frac{1}{r^6} + o\left(\frac{1}{r^6}\right), \qquad (31) $$

where the angular distribution is given by the functions $\gamma_{|m_e|}$ which read for the different values of $m_e$

$$ \gamma_0(\vartheta) = \sin^2 \vartheta, \qquad \gamma_1(\vartheta) = \frac{1}{4}(1 + \cos^2 \vartheta), \qquad (32) $$

and are represented in FIG. 1. We point out that the angular distribution we have found matches that of a classical dipole that is either oscillating linearly for $\gamma_0$ or circularly for $\gamma_1$ [42, 43]. We further remark that in the $m_e = 0$ case, the leading term of the asymptotics is zero for $\vartheta = 0$ due to the shape of the angular distribution. However, we emphasize that the spontaneous emission in this direction is not zero but features a faster falloff.

The result obtained in equation (31) exhibits an algebraic decrease for the energy density which is therefore far from the quasiexponential limit introduced in [16], and cannot be attributed only to the blurring of the light cone as described in [37]. This behavior is an illustration of the intrinsic nonlocality of the energy density of single-photon states, as shown for particular cases in [17] and for all single-photon states in [18]. It is also an other demonstration of the importance of multiphoton components that are predicted by the coupling to guarantee a causal emission [35].

## V. CONCLUSION AND OUTLOOK

Under the assumption that the spontaneous emission of a Hydrogen atom produces a single-photon, we have determined the spatial nonlocality properties, and shown that at any given time the expectation value of the local energy density is not zero anywhere, and decreases asymptotically like a power $1/r^6$ as a function of the distance $r$ of the detector to the atom. We remark that this power law decrease is far from the maximal decrease, which is
quasiexponential, as determined in [16]. It was shown in [18] that all single-photon states have a nonzero expectation value of the energy in any finite volume. Our result provides a concrete illustration of the intrinsic nonlocal property of single photons through a calculation that does not only consider the field but its source as well, and it shows that the nonlclity exists even before the photon state is completely emitted. From these results, one can conclude that the process of spontaneous emission cannot produce a strict single-photon state since it would contradict causality [19, 20]. It is known that there are states for which the expectation values of all observables are localized in space [11, 44–47]. These states contain necessarily n-photon components for arbitrarily high n, with weights that can be small but nonzero. Coherent states are an example. An open question that should be addressed is the determination of the weights of the n-photon components in concrete emission processes, like spontaneous emission or in the production of photons at demand by the different techniques that are being developed experimentally.

ACKNOWLEDGMENTS

We thank Jonas Lampart for fruitful discussions. This work was supported by the “Investissements d’Avenir” program, project ISITE-BFC/IQUINS (ANR-15-IDEX-03), QUACO-PRC (ANR-17-CE40-0007-01) and the EUR-EIPHI Graduate School (17-EURE-0002). We also acknowledge support from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 765075 (LIMQUET) and from France 2030 (ANR-22-CMAS-0001) QuanTEdu-France. M. Federico also acknowledges funding from European Union’s Horizon Europe research and innovation programme under the project Quantum Secure Network Partnership (QSNP, grant agreement No 101114043).

Appendix A: Position space representations and related operators

To analyze the spatial properties of photons, one needs to use a position space representation of the states. This
can be done using several equivalent choices [34]. In this work, we have chosen to work with the Landau-Peierls (LP)
representation which is constructed as follows

$$
\vec{\psi}(\vec{x}) = \sqrt{\frac{\epsilon_0}{2\hbar}} \left[ \Omega^{1/2} \vec{A}(\vec{x}) - i\Omega^{-1/2} \vec{E}(\vec{x}) \right], \quad (A1)
$$

where $\vec{A}$ is the classical potential vector in the Coulomb gauge and $\vec{E}$ is the classical electric field. The operators $\Omega^{\pm 1/2}$ are self-adjoint operators constructed from the frequency operator $\Omega = c(-\Delta)^{1/2}$ for the Laplacian $\Delta$. These operators are well defined and Refs [18, 34, 48] give the more details of their construction. The LP representation is made in the Hilbert space

$$
\mathcal{H}_{LP} = \{ \vec{\psi} \mid \nabla \cdot \vec{\psi} = 0, \langle \vec{\psi} | \vec{\psi}' \rangle_{LP} < \infty \}, \quad (A2)
$$

with the scalar product

$$
\langle \vec{\psi} | \vec{\psi}' \rangle_{LP} = \int_{\mathbb{R}^3} d^3 x \, \vec{\psi}^* \cdot \vec{\psi}', \quad (A3)
$$

and Maxwell’s equations written in this representation read

$$
i \frac{\partial \vec{\psi}}{\partial t} = \Omega \vec{\psi}, \qquad \nabla \cdot \vec{\psi} = 0. \tag{A4}
$$

This Hilbert space can then be used to construct the bosonic Fock space of states as

$$
\mathbb{F}^{\mathfrak{Q}}(\mathcal{H}_{LP}) = \bigoplus_{n=0}^{\infty} \mathcal{H}_{LP}^{\otimes s^n}, \qquad (A5)
$$

where $\mathcal{H}_{LP}^{\otimes s^n}$ is the symmetrized $n$-times tensor product of $\mathcal{H}_{LP}$ with itself. Annihilation and creation operators directly constructed on arbitrary pulse-shaped functions can then be used to write states of the fields from the vacuum $|\emptyset\rangle$:

$$
\hat{B}_{\vec{\psi}} |\emptyset\rangle = 0, \qquad \hat{B}_{\vec{\psi}}^{\dagger} |\emptyset\rangle = |\vec{\psi}\rangle. \tag{A6}
$$
They satisfy the general bosonic commutation relations [34, 49, 50]

$$[\hat{B}_{\vec{\psi}}, \hat{B}_{\vec{\psi}'}^{\dagger}] = \langle \vec{\psi} | \vec{\psi}' \rangle_{LP}. \qquad (A7)$$

We emphasize that the definition is done for any function of $\mathcal{H}_{LP}$. As usual, one can also use any generalized basis like the circularly polarized plane waves or the helicity vector spherical harmonics introduced in Appendix B, as we did in eq. (4). Defining a basis allows also to express the electromagnetic field operators. Indeed, if we consider a basis $\{\vec{\phi}_{\kappa}\}$, the classical LP function $\vec{\psi}(\vec{x})$ can be expanded as

$$\vec{\psi}(\vec{x}) = \sum_{\kappa} \vec{\phi}(\vec{x}) z_{\kappa}, \qquad (A8)$$

with $z_{\kappa} = \langle \vec{\phi}_{\kappa} | \vec{\psi} \rangle_{LP}$. Its quantum equivalent is thus

$$\vec{\Psi}(\vec{x}) = \sum_{\kappa} \vec{\phi}(\vec{x}) \hat{B}_{\vec{\phi}_{\kappa}}, \qquad (A9)$$

from which the real electromagnetic field observables are expressed as the following Hermitian operators

$$\vec{A}(\vec{x}) = \sqrt{\frac{\hbar}{2\epsilon_0}} \Omega^{-1/2} (\vec{\Psi}(\vec{x}) + \vec{\Psi}^{\dagger}(\vec{x})), \qquad \vec{E}(\vec{x}) = i\sqrt{\frac{\hbar}{2\epsilon_0}} \Omega^{1/2} (\vec{\Psi}(\vec{x}) - \vec{\Psi}^{\dagger}(\vec{x})). \qquad (A10)$$

Another important tool that is used is the helicity operator, defined as

$$\Lambda = c \Omega^{-1} \nabla \times . \qquad (A11)$$

For transverse fields, its spectrum is $\{\pm 1\}$ and one can thus decompose any transverse field $\vec{v}$ into a sum of a positive- and a negative-helicity part $\vec{v}^{(h\pm)}$,

$$\vec{v} = \vec{v}^{(h+)} + \vec{v}^{(h-)}, \qquad (A12)$$

where $\Lambda \vec{v}^{(h\pm)} = \pm \vec{v}^{(h\pm)}$. This splitting is an important property for the energy density mean value (21) since it leads to the nonlocality of the energy density of all single-photon states [18].

## Appendix B: Helicity vector spherical harmonics

To compute the state of the emitted photon, we have used the basis of helicity vector spherical harmonics $|k, J, M, \lambda\rangle \equiv \vec{\psi}_{k, J, M}^{(\lambda)}(\vec{x})$ that are eigenfunctions of the following operators [5, 51, 52]

$$-\Delta \vec{\psi}_{k, J, M}^{(\lambda)} = k^2 \vec{\psi}_{k, J, M}^{(\lambda)}, \qquad k \in (0, \infty), \qquad (B1a)$$

$$\hat{J}^2 \vec{\psi}_{k, J, M}^{(\lambda)} = J(J+1) \vec{\psi}_{k, J, M}^{(\lambda)}, \qquad J \in \{0, 1, 2, \dots\}, \qquad (B1b)$$

$$\hat{J}_3 \vec{\psi}_{k, J, M}^{(\lambda)} = M \vec{\psi}_{k, J, M}^{(\lambda)}, \qquad M \in \{-J, \dots, J\}, \qquad (B1c)$$

$$\Lambda \vec{\psi}_{k, J, M}^{(\lambda)} = \lambda \vec{\psi}_{k, J, M}^{(\lambda)}, \qquad \lambda = \pm 1, \qquad (B1d)$$

where

$$\vec{J} = \vec{L} + \vec{S}, \qquad (B2a)$$

$$\vec{L} = \nabla \times \vec{x} = \vec{x} \times \nabla, \qquad (B2b)$$

$$\vec{S} = (\hat{S}_1, \hat{S}_2, \hat{S}_3)^T, \qquad [\hat{S}_l]_{m,n} = -i\epsilon_{lmn}, \qquad (B2c)$$

where $l, m, n \in \{1, 2, 3\}$. The helicity vector spherical eigenfunctions can be written in spherical coordinates $(r, \vartheta, \varphi)$ as

$$\vec{\psi}_{k, J, M}^{(\lambda)}(r, \vartheta, \varphi) = \frac{i}{\sqrt{2}} \left[ \sqrt{\frac{J+1}{2J+1}} \vec{\psi}_{k, J, M}^{J-1}(r, \vartheta, \varphi) - \sqrt{\frac{J}{2J+1}} \vec{\psi}_{k, J, M}^{J+1}(r, \vartheta, \varphi) - i\lambda \vec{\psi}_{k, J, M}^{J}(r, \vartheta, \varphi) \right], \qquad (B3)$$
where the functions $\vec{\psi}_{k,J,M}^L$ can be expressed in terms of vector spherical harmonics $\vec{Y}_{J,M}^L$, i.e., eigenfunctions of $\hat{J}^2$, $\hat{J}_3$, $\hat{L}^2$ and $\hat{S}^2$ but not of $\Lambda$, and spherical Bessel functions $j_L(kr)$ as

$$ \vec{\psi}_{k,J,M}^L(r, \vartheta, \varphi) = \sqrt{\frac{2}{\pi}} k \, j_L(kr) \vec{Y}_{J,M}^L(\vartheta, \varphi). \qquad (B4) $$

We remark that by constructing the functions $\vec{\psi}_{k,J,M}^{(\lambda)}$ from $\vec{\psi}_{k,J,M}^L$ we obtain eigenfunctions of $-\Delta$, $\hat{J}^2$, $\hat{J}_3$ and $\Lambda$ but not of $\hat{L}^2$ while $\vec{\psi}_{k,J,M}^L$ is an eigenfunction of $-\Delta$, $\hat{J}^2$, $\hat{J}_3$ and $\hat{L}^2$ but not of $\Lambda$.

The helicity spherical vector eigenfunctions $\{\vec{\psi}_{k,J,M}^{(\lambda)}\}$ are transverse, $\nabla \cdot \vec{\psi}_{k,J,M}^{(\lambda)} = 0$, orthonormal

$$ \int_0^\infty dr \, r^2 \int_0^\pi d\vartheta \sin\vartheta \int_0^{2\pi} d\varphi \, \vec{\psi}_{k,J,M}^{\lambda*} \cdot \vec{\psi}_{k',J',M'}^{\lambda'} = \delta(k-k')\delta_{J,J'}\delta_{M,M'}\delta_{\lambda,\lambda'}, \qquad (B5) $$

and form a complete set of the subspace of transverse field since the set of function $\{\vec{\psi}_{k,J,M}^{J-1}, \vec{\psi}_{k,J,M}^J, \vec{\psi}_{k,J,M}^{J+1}\}$ from which they are defined is a complete set with

$$ \int_0^\infty d^3k \sum_{J,M,L} \left[ \vec{\psi}_{k,J,M}^L(r', \vartheta', \varphi') \right]_L \left[ \vec{\psi}_{k,J,M}^{L*}(r, \vartheta, \varphi) \right]_n = \frac{1}{r^2 \sin\vartheta} \delta_{l,n} \delta(r-r') \delta(\vartheta-\vartheta') \delta(\varphi-\varphi'). \qquad (B6) $$

The last equation is a consequence of the completeness of the Bessel functions $j_L(kr)$ and vector spherical harmonics $\vec{Y}_{J,M}^L$ [51].

## Appendix C: Asymptotic analysis of radial integrals

In this appendix, we compute the asymptotic behavior of the radial integrals that appear in the calculation of the energy density of a spontaneously emitted photon from the Lyman-α transition of a Hydrogen atom. The integrals have the following forms

$$ F_0(r,t) = \left(\frac{2}{3}\right)^3 \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q^2}{(1+q^2)^2} j_0(qr'), \qquad (C1a) $$

$$ F_1(r,t) = -i\lambda \left(\frac{2}{3}\right)^{5/2} \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q^2}{(1+q^2)^2} j_1(qr'), \qquad (C1b) $$

$$ F_2(r,t) = - \left(\frac{2}{3}\right)^3 \sqrt{\frac{\alpha^5}{2c}} \frac{mc^2}{\pi \hbar r_B} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q^2}{(1+q^2)^2} \left( \frac{3}{qr'} j_1(qr') - j_0(qr') \right), \qquad (C1c) $$

where we have introduced the dimensionless parameters

$$ A = \frac{\tilde{\Gamma}_a}{cK} = \frac{\Gamma_a}{2cK}, \qquad B = \frac{\Omega_a}{cK} = \frac{\omega_a + \Delta_a}{cK}, \qquad (C2a) $$

$$ q = \frac{k}{K}, \qquad p = cKt, \qquad r' = Kr. \qquad (C2b) $$

To obtain the asymptotic behavior of the $F_L$ integrals, and due to the expression of the spherical Bessel functions

$$ j_0(qr') = \frac{\sin(qr')}{qr'}, \qquad j_1(qr') = \frac{\sin(qr')}{q^2 r'^2} - \frac{\cos(qr')}{qr'}, \qquad (C3) $$

one needs to analyze integrals of the form

$$ \int_0^\infty dq R(q) \sin(qr'), \qquad \int_0^\infty dq R(q) \cos(qr'). \qquad (C4) $$
Here, $R(q)$ is a function with the properties that for some $N_s \ge 0$ even or $N_c \ge 0$ odd,

$$ \left. \frac{d^n R(q)}{dq^n} \right|_{q=0} = \begin{cases} 0 & \text{for } n < N_{c,s}, n \text{ same parity as } N_{c,s}, \\ \text{cst} \neq 0 & \text{for } n = N_{c,s}, n \text{ same parity as } N_{c,s}, \end{cases} \qquad (C5a) $$

$$ \lim_{q \to \infty} \frac{d^n R(q)}{dq^n} = 0, \qquad (C5b) $$

and

$$ \int_{0}^{\infty} dq \left| \frac{d^n R(q)}{dq^n} \right| < \infty, \quad \text{for all } n < N_{c,s}. \qquad (C6) $$

The indices $c, s$, are needed to treat the cosine and sine integrals in (C4), respectively. Using these properties and some integrations by parts the sine integral becomes

$$ \int_{0}^{\infty} dq R(q) \sin(qr') = -\frac{1}{r'} \left[ R(q) \cos(qr') \right]_{0}^{\infty} + \frac{1}{r'} \int_{0}^{\infty} dq \frac{dR(q)}{dq} \cos(qr') \qquad (C7a) $$

$$ = \frac{1}{r'} R(0) - \frac{1}{r'} \int_{0}^{\infty} dq \frac{dR(q)}{dq} \cos(qr') \qquad (C7b) $$

$$ = \frac{1}{r'} R(0) + o \left( \frac{1}{r'} \right). \qquad (C7c) $$

The last equality is obtained by the Riemann-Lebesgue lemma [53] which states that

$$ \left| \int_{0}^{\infty} dq \frac{dR(q)}{dq} \cos(qr') \right| \underset{r' \to +\infty}{\longrightarrow} 0, \qquad (C8) $$

if $\int_0^\infty dq \left|\frac{dR(q)}{dq}\right| < \infty$, which is true according to (C6). Therefore, the first term of (C7c) gives the asymptotics for large $r'$ if $N_s = 0$, otherwise, one has to iterate the integration by parts as

$$ \int_{0}^{\infty} dq R(q) \sin(qr') = \frac{1}{r'^2} \left[ \frac{dR}{dq} \sin(qr') \right]_{0}^{\infty} - \frac{1}{r'^2} \int_{0}^{\infty} dq \frac{d^2 R(q)}{dq^2} \sin(qr') \qquad (C9a) $$

$$ = -\frac{1}{r'^2} \int_{0}^{\infty} dq \frac{d^2 R(q)}{dq^2} \sin(qr') \qquad (C9b) $$

$$ = \frac{1}{r'^3} \left[ \frac{d^2 R}{dq^2} \cos(qr') \right]_{0}^{\infty} + \frac{1}{r'^3} \int_{0}^{\infty} dq \frac{d^3 R(q)}{dq^3} \cos(qr') \qquad (C9c) $$

$$ = \frac{1}{r'^3} \left. \frac{d^2 R}{dq^2} \right|_{q=0} + \frac{1}{r'^3} \int_{0}^{\infty} dq \frac{d^3 R(q)}{dq^3} \cos(qr') \qquad (C9d) $$

$$ = \frac{1}{r'^3} \left. \frac{d^2 R}{dq^2} \right|_{q=0} + o \left( \frac{1}{r'^3} \right), \qquad (C9e) $$

and again, the first term gives the asymptotics by the Riemann-Lebesgue lemma if $N_s = 2$. Otherwise, the same process can be continued until one reaches $N_s$. A similar process can be done for the cosine integral and gives after one step

$$ \int_{0}^{\infty} dq R(q) \cos(qr') = -\frac{1}{r'^2} \left. \frac{dR}{dq} \right|_{q=0} + \frac{1}{r'^2} \int_{0}^{\infty} dq \frac{d^2 R(q)}{dq^2} \cos(qr') \qquad (C10a) $$

$$ = -\frac{1}{r'^2} \left. \frac{dR}{dq} \right|_{q=0} + o \left( \frac{1}{r'^2} \right), \qquad (C10b) $$

by the Riemann-Lebesgue lemma. This step is enough if $N_c = 1$ but others can be performed if not.

To apply this technique, we now need the concrete form of the functions $R(q)$ and thus we rewrite the $F_L$ integrals as

$$ F_0 = \left(\frac{2}{3}\right)^3 \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \frac{1}{r'} \int_0^\infty dq \left(R_+^{(1)}(q) + R_-^{(1)}(q)\right) \sin(qr'), \qquad (C11a) $$

$$ F_1 = -i\lambda \left(\frac{2}{3}\right)^{5/2} \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \left\{ \frac{1}{r'^2} \int_0^\infty dq \left(R_+^{(0)}(q) + R_-^{(0)}(q)\right) \sin(qr') - \frac{1}{r'} \int_0^\infty dq \left(R_+^{(1)}(q) + R_-^{(1)}(q)\right) \cos(qr') \right\}, \qquad (C11b) $$
where the R functions are

$$R_{+}^{(1)}(q) = \frac{e^{-iqp}}{A + i(B - q)} \frac{q}{(1 + q^2)^2}, \qquad R_{-}^{(1)}(q) = -\frac{e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q}{(1 + q^2)^2}, \qquad (\text{C12a})$$

$$R_{+}^{(0)}(q) = \frac{e^{-iqp}}{A + i(B - q)} \frac{1}{(1 + q^2)^2}, \qquad R_{-}^{(0)}(q) = -\frac{e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{1}{(1 + q^2)^2}. \qquad (\text{C12b})$$

The integral $F_2$ can be deduced rapidly from the other two as we will discuss later. The asymptotics for $F_0$ and $F_1$ can be computed using the iterative process and we obtain

$$F_0(r,t) \underset{r' \to \infty}{\sim} \left(\frac{2}{3}\right)^3 \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \frac{1}{r'^4} \left( \left. \frac{d^2 R_+^{(1)}}{dq^2} \right|_{q=0} + \left. \frac{d^2 R_-^{(1)}}{dq^2} \right|_{q=0} \right), \qquad (\text{C13a})$$

$$F_1(r,t) \underset{r' \to \infty}{\sim} -i\lambda \left(\frac{2}{3}\right)^{5/2} \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \frac{1}{r'^3} \left( R_+^{(0)}(0) + R_-^{(0)}(0) + \left. \frac{d^2 R_+^{(1)}}{dq^2} \right|_{q=0} + \left. \frac{d^2 R_-^{(1)}}{dq^2} \right|_{q=0} \right), \qquad (\text{C13b})$$

which can be further simplified using

$$R_{+}^{(0)}(q=0) = \frac{1}{A + iB}, \qquad R_{-}^{(0)}(q=0) = -\frac{e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + iB}, \qquad (\text{C14a})$$

$$\left. \frac{d^2 R_{+}^{(1)}}{dq^2} \right|_{q=0} = \frac{i}{(A + iB)^2} - \frac{ip}{A + iB}, \qquad \left. \frac{d^2 R_{-}^{(1)}}{dq^2} \right|_{q=0} = -i \frac{e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{(A + iB)^2}. \qquad (\text{C14b})$$

Since $F_0 \sim F_1/r'$, we compute only $F_1$ which will dominate for large $r'$. We obtain then the result used in Section IV by inserting (C14) in (C13b)

$$F_1(r', t) \underset{r' \to \infty}{\sim} -i\lambda \left(\frac{2}{3}\right)^{5/2} \sqrt{\frac{\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \mathcal{T}(t) \frac{1}{r'^3}, \qquad (\text{C15})$$

where the dimensionless function of time reads

$$\mathcal{T}(t) = \frac{1}{A + iB} \left( 1 - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t} - 2icKt + i \frac{1 - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + iB} \right). \qquad (\text{C16})$$

Regarding $F_2$ we recall its expression

$$F_2(r,t) = -\left(\frac{2}{3}\right)^3 \sqrt{\frac{\alpha^5}{2c}} \frac{mc^2}{\pi \hbar r_B} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q^2}{(1 + q^2)^2} \left( \frac{3}{qr'} j_1(qr') - j_0(qr') \right), \qquad (\text{C17})$$

from which we see that the term involving $j_0$ will behave as $F_0$ and can thus be neglected. For the term involving $j_1$, one has to be more cautious since it does not have the same form as what we had for $F_1$. It yields

$$F_2(r', t) \underset{r' \to \infty}{\sim} -\left(\frac{2}{3}\right)^2 \sqrt{\frac{2\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \frac{1}{r'} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{q}{(1 + q^2)^2} j_1(qr') \qquad (\text{C18a})$$

$$= -\left(\frac{2}{3}\right)^2 \sqrt{\frac{2\alpha^5}{c}} \frac{mc^2}{\pi \hbar r_B} \left\{ \frac{1}{r'^3} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{1}{q(1 + q^2)^2} \sin(qr') \right. \\ \left. - \frac{1}{r'^2} \int_0^\infty dq \frac{e^{-iqp} - e^{-i(\Omega_a - i\tilde{\Gamma}_a)t}}{A + i(B - q)} \frac{1}{(1 + q^2)^2} \cos(qr') \right\}, \qquad (\text{C18b})$$

from which one can already see that $F_2$ will decrease faster than $F_1$. Indeed, the first term scales like $1/r'^3$ before making the iterative process to analyze the integral. Moreover, the function of $q$ to integrate does not yield a constant when $q \to 0$ meaning that the iterative process must be done more than one time. Thus, this term will decrease faster than $1/r'^3$ and can be neglected. The second term involves a cosine for which the first iterative step (C10) already
scales as $1/r'^2$, so that the full term would scale like $1/r'^4$, decreasing faster than $F_1$. The leading term that is needed to compute the mean value is thus given by $F_1$ and if one reinstates the units using (C2), it yields

$$F_1(r,t) \sim -i\lambda \left(\frac{2}{3}\right)^{11/2} \sqrt{\frac{\alpha^5}{c} \frac{mc^2 r_B^2}{\pi\hbar}} \mathcal{T}(t) \frac{1}{r^3}. \quad (C19)$$

[1] V. Weisskopf and E. Wigner, Zeitschrift für Physik **63**, 54 (1930).

[2] H. Spohn, *Dynamics of Charged Particles and Their Radiation Field* (Cambridge University Press, 2009).

[3] L. Novotny and B. Hecht, *Principles of Nano-Optics* (Cambridge University Press, 2012).

[4] J. Dereziński and R. Früboes, Reports on Mathematical Physics **50**, 433 (2002).

[5] H. E. Moses, Physical Review A **8**, 1710 (1973).

[6] H. E. Moses, Physical Review A **11**, 2196 (1975).

[7] H. E. Moses, Physical Review A **11**, 387 (1975).

[8] E. Grimm and V. Ernst, Journal of Physics A: Mathematical, Nuclear and General **7**, 1664 (1974).

[9] E. Grimm and V. Ernst, Zeitschrift für Physik A Atoms and Nuclei **274**, 293 (1975).

[10] J. Seke and W. N. Herfort, Physical Review A **38**, 833 (1988).

[11] J. M. Knight, Journal of Mathematical Physics **2**, 459 (1961).

[12] G. C. Hegerfeldt, Physical Review D **10**, 3320 (1974).

[13] G. C. Hegerfeldt, Physical Review Letters **72**, 596 (1994).

[14] G. C. Hegerfeldt, arXiv:quant-ph/9707016 (1997).

[15] G. C. Hegerfeldt, in *Irreversibility and Causality Semigroups and Rigged Hilbert Spaces* (Springer Berlin Heidelberg, 1998) pp. 238–245.

[16] I. Bialynicki-Birula, Physical Review Letters **80**, 5247 (1998).

[17] I. Bialynicki-Birula and Z. Bialynicka-Birula, Physical Review A **79**, 032112 (2009).

[18] M. Federico and H. R. Jauslin, Physical Review A **108**, 043720 (2023).

[19] P. R. Berman, Physical Review A **69**, 022101 (2004).

[20] I. Dolce, R. Passante, and F. Persico, Physics Letters A **355**, 152 (2006).

[21] J. Gulla and J. Skaar, Physical Review Letters **126**, 073601 (2021).

[22] K. Ryen, J. Gulla, and J. Skaar, International Journal of Theoretical Physics **61**, 10.1007/s10773-022-05133-7 (2022).

[23] J. Gulla, K. Ryen, and J. Skaar, Physical Review A **108**, 063708 (2023).

[24] A. K. Biswas, G. Compagno, G. M. Palma, R. Passante, and F. Persico, Physical Review A **42**, 4291 (1990).

[25] D. Buchholz and J. Yngvason, Physical Review Letters **73**, 613 (1994).

[26] O. Keller, Physical Review A **62**, 022111 (2000).

[27] K. W. Chan, C. K. Law, and J. H. Eberly, Physical Review Letters **88**, 100402 (2002).

[28] M. V. Fedorov, M. A. Efremov, A. E. Kazakov, K. W. Chan, C. K. Law, and J. H. Eberly, Physical Review A **72**, 032110 (2005).

[29] V. Debierre, T. Durt, A. Nicolet, and F. Zolla, Physics Letters A **379**, 2577 (2015).

[30] V. Debierre, I. Goessens, E. Brainis, and T. Durt, Physical Review A **92**, 023825 (2015).

[31] V. Debierre and T. Durt, Physical Review A **93**, 023847 (2016).

[32] V. Debierre, B. Stout, and T. Durt, arXiv:1809.11070 (2018).

[33] I. Bialynicki-Birula, in *Progress in Optics*, edited by E. Wolf (Elsevier, 1996).

[34] M. Federico and H. R. Jauslin, Journal of Physics A: Mathematical and Theoretical **56**, 235302 (2023).

[35] G. Compagno, R. Passante, and F. Persico, Journal of Modern Optics **37**, 1377 (1990).

[36] G. Compagno, R. Passante, and F. Persico, *Atom-Field Interactions and Dressed Atoms* (Cambridge University Press, 1995).

[37] M. I. Shirokov, Soviet Physics Uspekhi **21**, 345 (1978).

[38] L. Landau and R. Peierls, Zeitschrift für Physik **62**, 188 (1930).

[39] L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, 1995).

[40] K. O. Friedrichs, Communications on Pure and Applied Mathematics **1**, 361 (1948).

[41] T. D. Lee, Physical Review **95**, 1329 (1954).

[42] J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (John Wiley & Sons Inc, 1999).

[43] D. A. Steck, *Quantum and Atom Optics* (Available online at http://steck.us/teaching, 2023) revision 0.13.21.

[44] A. L. Licht, Journal of Mathematical Physics **4**, 1443 (1963).

[45] A. L. Licht, Journal of Mathematical Physics **7**, 1656 (1966).

[46] S. De Bièvre, in *Large Coulomb Systems*, Springer Berlin Heidelberg, 2006) pp. 15–61.

[47] S. De Bièvre, in *Contributions in Mathematical Physics* (Hindustan Book Agency, 2007) pp. 123–146.

[48] M. Federico, V. Dorier, S. Guérin, and H. R. Jauslin, Journal of Physics B: Atomic, Molecular and Optical Physics **55**, 174002 (2022).

[49] J. C. Garrison and R. Y. Chiao, *Quantum Optics* (Oxford University Press, 2008).

[50] C. Fabre and N. Treps, Reviews of Modern Physics **92**, 035005 (2020).
[51] D. A. Varshalovich, *Quantum Theory of Angular Momentum* (World Scientific Pub., 1988).

[52] L. Dai, M. Kamionkowski, and D. Jeong, Physical Review D **86**, 125013 (2012).

[53] W. Appel, *Mathematics for Physics and Physicists* (Princeton University Press, 2007) p. 672.