**Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A. Classical Antenna Theory**

Let the source or current Green's function be the dyadic tensor **<sup>F</sup>**(**<sup>x</sup>**, **<sup>x</sup>**, *t* − *t*) [28,46–48,89]. The forward Green's function is the standard retarded radiation Green's function **<sup>G</sup>**(**x** − **<sup>x</sup>**, *t* − *t*) [49]. Then, we have the following two fundamental relations [28,90]:

$$\mathbf{J}(\mathbf{x},t) = \int\_{S} \mathbf{d}^2 \mathbf{x}' \, \mathbf{d}t' \, \mathbb{E}(\mathbf{x}, \mathbf{x}', t - t') \cdot \mathbf{E}^{\infty}(\mathbf{x}', t'), \tag{A1}$$

$$\mathbf{E}(\mathbf{x},t) = \int\_{S} \mathbf{d}^3 \mathbf{x}' \, \mathrm{d}t' \, \mathbf{\overline{G}}(\mathbf{x}, \mathbf{x}', t - t') \cdot \mathbf{J}(\mathbf{x}', t'), \tag{A2}$$

where **<sup>E</sup>**ex(**<sup>x</sup>**, *t*) is the excitation electric field, **J**(**<sup>x</sup>**, *t*) is the radiating antenna current on its compact surface *S* supporting the electromagnetic boundary conditions, and **<sup>E</sup>**(**<sup>x</sup>**, *t*) is the radiated field. For simplicity, we assume perfect-electric conducting (PEC) antennas where the magnetic field does not contribute to the current Green's function. Another radiating Green's functions similar to **<sup>G</sup>**(**<sup>x</sup>**, **<sup>x</sup>**, *t* − *t*) is needed in expressions such as (A2) in order to obtain the magnetic field **<sup>B</sup>**(**<sup>x</sup>**, *t*), which is always present in any radiation problem beside the electric field; see [89–91] for more details and applications.

The relation (A1) captures the antenna's fundamental Mode A, where a source field **<sup>E</sup>**ex(**<sup>x</sup>**, *t*) produces a radiating current via the current Green's function **<sup>F</sup>**(**<sup>x</sup>**, **<sup>x</sup>**, *t* − *<sup>t</sup>*), which in turn generates the radiated fields **<sup>E</sup>**(**<sup>x</sup>**, *t*) and **<sup>B</sup>**(**<sup>x</sup>**, *t*) throughout the region exterior to the surface *S*. To complete the antenna-based wireless communication system, a third Mode C, where a receiving (Rx) antenna, placed at some distance from *S*, interacts with the radiated field in order to produce an observable receive port current **J**rx(**x**) by means of the formula

$$\mathbf{J}\_{\rm rx}(\mathbf{x},t) = \int\_{S\_{\rm rx}} \mathbf{d}^2 \mathbf{x}' \mathbf{d}t' \,\overline{\mathbf{F}}\_{\rm rx}(\mathbf{x}, \mathbf{x}', t - t') \cdot \mathbf{E}(\mathbf{x}', t'), \tag{A3}$$

where **<sup>F</sup>**rx(**<sup>x</sup>**, **<sup>x</sup>**, *t* − *t*) is the receive current Green's function, which is generally different from the transmit (Tx) current Green's function **<sup>F</sup>**(**<sup>x</sup>**, **<sup>x</sup>**, *t* − *<sup>t</sup>*). The three relations and Green's functions in (A1)–(A3) fully describe the antenna system in classical electromagnetic theory [92], with obvious generalization to magnetic field interactions included in essentially the same logic [90].

### **Appendix B. The Relativistic Four-Vector Formalism**

In special relativity, everything takes place in the four-dimensional Minkowski spacetime, denoted by M4, a linear vector space with a special metric, the Lorentz metric, given by

$$\mathcal{g}\_{\mu\nu} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix} . \tag{A4}$$

We follow the standard convention of denoting spacetime indices using Greek letters *μ* = 0, 1, 2, 3, with the time-like component always given the 0th index and the spacelike part denoted by Roman letter. Four vectors with upper and lower indices are dual vectors. The converse of the metric tensor *gμν* is *gμν*, i.e., *<sup>g</sup>μνgμν* = *δνμ*, where *δνμ* is the Kronecker delta function. Note that summation over repeated indices is implied (the Einstein repeated index convention). In general, we raise or lower indices by contracting with the fundamental tensor in relations such as *<sup>a</sup>μ* = *<sup>g</sup>μνa<sup>ν</sup>* or *aμ* = *gμνaν*. From this we can see, for example, that *a*0 = *a*0, and *ai* = <sup>−</sup>*ai*. The dot product between two four vectors *aμ* and *bμ* is defined as *<sup>a</sup>μbμ* := *<sup>a</sup>μgμνb<sup>ν</sup>* = *<sup>g</sup>μνaμb<sup>ν</sup>*. We also write this as *a* · *b* := *<sup>a</sup>μbμ* where *a* and *b* stands for the four-vectors *aμ* and *bμ*. In particular, the all-important inner product *p* · *x* between four-momentum *pμ* and spacetime position *<sup>x</sup>μ*, given by

$$p \cdot \mathbf{x} = p^{\mu} \mathbf{x}\_{\mu} = \omega t - \mathbf{k} \cdot \mathbf{x}\_{\prime} \tag{A5}$$

will be used frequently. Note that in M4, the proper "differential operators four-vectors" are given by *∂μ* = (*∂*/*∂t*, −∇) and *∂μ* = (*∂*/*∂t*, ∇). Table A1 provides a compact summary of the main relativistic formulas.



### **Appendix C. Natural Units**

In relativistic quantum field theory, it is customary to use a system of units in which *h*¯ = 1, *c* = 1. Here, the action *S* is dimensionless. Mass, energy, wavevector, and momentum have the same dimension, which is the inverse of the dimension of length and time:

$$[\text{mass}] = [\text{energy}] = [\text{momentum}] = [\text{length}]^{-1} = [\text{time}]^{-1} \tag{A6}$$

Therefore, the physics of the world is reduced to measurements in two units, length and time. Moreover, in four dimensions, the creation and annihilation operators have the dimension of an inverse energy, i.e., [*<sup>a</sup>***p**]=[*a*†**p**]=[energy]−1.

### **Appendix D. Dirac Interaction Picture**

Consider a Hamiltonian *H* that can be written as

$$H = H\_0 + H'(t),\tag{A7}$$

where *H*0 is the *free* Hamiltonian of the system (time-independent function), while *H*(*t*) is the *interaction* (time-dependent) Hamiltonian, which is usually assumed to be turned on at a specific time point. In the Dirac picture, our goal is to *decouple* the free dynamics of the system, i.e., that governed by a Hamiltonian of the form *H* = *H*0 (no interaction or *H* = 0), from the pure interaction component, i.e., that which is solely due to the term *H* in the actual full Hamiltonian (A7). Now, let us define the following time-dependent operator

$$O\_{\rm I}(t) := e^{iH\_0t} O\_{\rm S} e^{-iH\_0t} \,. \tag{A8}$$

where *O*S is the Schrodinger picture (time-independent) operator. Differentiating (A8) with respect to time, we obtain the Heisenberg equation of motion

$$\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}O\_{\mathrm{I}}(t) = [O\_{\prime}H\_{0}].\tag{A9}$$

Therefore, the interaction picture operator *<sup>O</sup>*I(*t*) evolves in time according to the exact same law corresponding to the case when the full Hamiltonian is just the free part, i.e., *H* = *H*0.

What about the evolution of the interaction picture state, which will be denoted by *ψ*I(*t*)? To find out, note that all different pictures, the Schrodinger, Heisenberg, and Dirac, must agree on probability amplitudes. Therefore, for any two pairs of Schrodinger and Dirac states *ψ*S(*t*), *ϕ*S(*t*), and *ψ*I(*t*), *ϕ*I(*t*), respectively, the matrix elements in both pictures must be identical. Formally, we express this by the condition

$$
\langle \psi\_{\sf S}(t) | O\_{\sf S} \varphi\_{\sf S}(t) \rangle = \langle \psi\_{\sf I}(t) | O\_{\sf I}(t) \varphi\_{\sf I}(t) \rangle \\
= \left\langle \psi\_{\sf I}(t) e^{-iH\_0 t} \Big| O\_{\sf S} e^{-iH\_0 t} \varphi\_{\sf I}(t) \right\rangle. \tag{A10}
$$

Clearly, this implies that we should define the Dirac (interaction) state by

$$|\psi\_{\mathbf{l}}(t)\rangle := e^{\mathbf{l}^{\mathrm{H}\_{0}\mathbf{l}}}|\psi\_{\mathbf{\tilde{s}}}(t)\rangle. \tag{A11}$$

To find out how this state evolves, we differentiate with respect to time, obtaining

$$\frac{\mathbf{d}}{\mathbf{d}t} |\psi\_{\mathbf{l}}(t)\rangle = \left(-H\_{\mathbf{0}} + \mathbf{i}\frac{\mathbf{d}}{\mathbf{d}t}\right) |\psi\_{\mathbf{S}}(t)\rangle = H' |\psi\_{\mathbf{S}}(t)\rangle \tag{A12}$$

where we have made use of the Schrodinger equation

i

$$\operatorname{i}\frac{\mathbf{d}}{\mathbf{d}t}|\psi\_{\mathbf{S}}(t)\rangle = H|\psi\_{\mathbf{S}}(t)\rangle\tag{A13}$$

and the form (A7). Finally, using the definition (A11) of the interacting picture state, the relation (A12) becomes

$$\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}|\psi\_{\mathrm{I}}(t)\rangle = H\_{\mathrm{I}}(t)|\psi\_{\mathrm{I}}(t)\rangle,\tag{A14}$$

where the Dirac picture interacting Hamiltonian *H*I is given by

$$H\_{\rm I}(t) := e^{iH\_0t} H'(t) e^{-iH\_0t}. \tag{A15}$$

Let us summarize the main features of this Dirac picture:


3. In the Dirac picture, both state and operators evolve with time. However, the time evolution is decoupled into two distinct and independent components. First, all interaction (Dirac) picture operators evolve according to the free Hamiltonian as per the corresponding Heisenberg Equation (A9). Second, the Dirac state |*ψ*I(*t*) evolves independently according to the dynamic Equation (A14).

Therefore, the main advantage of the Dirac picture is that we can concentrate on the essential aspects of interaction as encoded in the Schrodinger-like dynamic Equation (A14). In this equation, Schrodinger or Heisenberg states and operators can be obtained using the usual transformation formula. On the other hand, note that operators in the Dirac picture, including in particular the interaction Hamiltonian (A15), all evolve according to the free evolution law (A9). In this case, this becomes

$$\mathbf{i}\frac{\mathbf{d}}{\mathbf{d}t}H\_{\mathbf{l}}(t) = [H\_{\mathbf{l}}(t), H\_{\mathbf{0}}]\_{\mathbf{l}} \tag{A16}$$

which says that all quantum fields inside *H* are to be evolved under the free Hamiltonian *H*0 in order to obtain the dynamics of *<sup>H</sup>*I(*t*).

Let us finally solve the dynamical evolution Equation (A14). To do so, introduce the *time-ordering symbol* T defined by

$$\mathcal{T}[A(t\_1)B(t\_2)] := \begin{cases} A(t\_1)B(t\_2), & t\_1 > t\_{2\prime} \\ B(t\_2)A(t\_1), & t\_2 > t\_1. \end{cases} \tag{A17}$$

Then, the general solution of (A14) can be expressed in terms of the *evolution operator <sup>P</sup>*(*<sup>t</sup>*2, *<sup>t</sup>*1), which evolves that initial state |*ψ*(*<sup>t</sup>*1) to the final state |*ψ*(*<sup>t</sup>*2) via the operator relation

$$|\psi\_{\mathbf{l}}(t\_2)\rangle = P(t\_2, t\_1)|\psi\_{\mathbf{l}}(t\_1)\rangle,\tag{A18}$$

where

$$P(t\_2, t\_1) := \mathcal{T} e^{-\mathrm{i} \int\_{t\_1}^{t\_2} \mathrm{d}t H\_1(t)}.\tag{A19}$$

The reader may verify that (A18) and (A19) indeed solves (A14) by direct substitution. Note that all operators inside the ordered exponential symbol T commute with each other.

### **Appendix E. On the Background to Theorem**

The complete and most general proof of Theorem 1 can be obtained by utilizing a generalized framework, that of *algebraic quantum field theory* [37,78]. The expansion (15) can be shown to be derivable from a suitable perturbative algebraic quantum field theory, e.g., the recent approach [93]. The advantage of choosing such a method is that one does not need to assume a concrete Lagrangian from the beginning but rather proceed to work directly with algebras of quantum fields and then use the structure of these algebras in order to construct the entire theory, including the quantum states themselves, which are generated internally. However, since the mathematical details are extensive, the full treatment will be presented in a separate paper. Nevertheless, note that special cases of Theorem 1 have already appeared repeatedly in the literature, though in quite different applications and within distinct contexts. For example, a special case of (15) seems to have been discussed by Coleman in his analysis of perturbation calculations in scattering theory [38]. Moreover, a less general form of the expansion (15), known as the Volterra series, is often presented in several textbooks on QFT when discussing the evaluation of the vacuum expectation persistence function in terms of the higher-point Green's (correlation) functions of the quantum field, e.g., see [38,54–58,67].

### **Appendix F. The Neutral Klein–Gordon Field Theory**

A massive neutral Klein–Gordon field with mass *m* can be fully captured by a real scalar field *φ*(*x*) whose Lagrangian density is given by [38,54,56,75]

$$\mathcal{L}(\mathbf{x}) = \frac{1}{2} \partial^{\mu} \partial\_{\mu} \phi(\mathbf{x}) + \frac{1}{2} m^{2} \phi(\mathbf{x})^{2},\tag{A20}$$

with the corresponding action integral being

$$S = \int \mathbf{d}^4 \mathbf{x} \,\mathcal{L}(\mathbf{x}).\tag{A21}$$

The Euler–Lagrange equation of motion

$$\frac{\delta S}{\delta \phi(\mathbf{x})} = 0 \tag{A22}$$

yields the field equation

$$
\partial^{\mu}\partial\_{\mu}\phi(\mathbf{x}) + m^2 = 0,\tag{A23}
$$

which is the relativistic wave equation of a particle the Klein–Gordon equation [38,75].

Note that within the framework of standard quantum field theory (QFT), fields are promoted to *operators* expanded in terms of plane-wave modes of the form

$$\text{operator} \times \exp\left(-\text{i}p\_{\mu}\mathbf{x}^{\mu}\right),\tag{A24}$$

where the "operator" is either creation *a*†**p** or annihilation *<sup>a</sup>***p** operator (for particles or antiparticles). Note further that by convention, the plane wave

$$\exp\left(-\mathrm{i}p\_{\mu}\mathbf{x}^{\mu}\right) = \exp\left[-\mathrm{i}(E\_{\mathbf{P}}t - \mathbf{p}\cdot\mathbf{x})\right] \tag{A25}$$

is taken to encode an *incoming* wave/particle with momentum **p**/wavevector **k** and energy *<sup>E</sup>***p**/frequency *ω* [38]. In the unified language of QFT, we say that the plane wave encodes a particle/antiparticle in a pure momentum state **p**. Within this convention, the energy/frequency of a particle/wave is *always positive*.<sup>25</sup>

By applying the canonical quantization algorithm [38] (see review in Appendix G), the quantum field *φ*(*x*) may be expanded into a continuous sum of spacetime modes (plane waves) as follows

$$\boldsymbol{\phi}(\mathbf{x}) = \int\_{\mathbf{P} \in \mathbb{R}^3} \frac{\mathbf{d}^3 \mathbf{p}}{(2\pi)^{3/2}} \frac{1}{(2\omega\_\mathbf{p})^{1/2}} \left[ a\_\mathbf{p} e^{-i p\_\mu \mathbf{x}^\mu} + a\_\mathbf{p}^\dagger e^{i p\_\mu \mathbf{x}^\mu} \right],\tag{A26}$$

with the dispersion relation (in natural units)

$$E\_{\mathbf{p}} = \omega\_{\mathbf{p}} = \omega\_{\mathbf{k}} = +\sqrt{|\mathbf{p}|^2 + m^2} = +\sqrt{|\mathbf{k}|^2 + m^2}.\tag{A27}$$

The creation and annihilation operators *a*†**p** and *<sup>a</sup>***p**, respectively, obey the standard canonical commutation relations

$$[a\_{\mathbf{p}'}a\_{\mathbf{p}'}] = 0,\ [a\_{\mathbf{p}'}^\dagger a\_{\mathbf{p}'}^\dagger] = 0,\ [a\_{\mathbf{p}'}a\_{\mathbf{p}'}^\dagger] = \delta(\mathbf{p} - \mathbf{p}'),\tag{A28}$$

where *δ* is the Dirac delta function.

### **Appendix G. The Relativistic Field-Theoretic Canonical Quantization Algorithm**

The general canonical quantization algorithm is shown in Algorithm A1, where we leave the nature of the field (scalar, vector, spin type, tensor, etc.) unspecified. In what follows, the detailed quantization algorithm applies to our main type of fields in this paper,

i.e., the massive spin-0 scalar field *φ*(*x*). Recall that in nonrelativistic quantum theory, a generic operator can be Fourier expanded in terms of creation and annihilation operators *a*†**p** and *<sup>a</sup>***p**, respectively. For instance, a generic scalar quantum field *φ*(*x*) is expected to be written as

$$\boldsymbol{\Phi}(\mathbf{x}) = \int\_{\mathbf{P} \in \mathbb{R}^3} \frac{\mathbf{d}^3 p}{(2\pi\mathbf{r})^3} \left[ a\_\mathbf{P} e^{-i p\_\mu \mathbf{x}^\mu} + a\_\mathbf{P}^\dagger e^{i p\_\mu \mathbf{x}^\mu} \right]. \tag{A29}$$

Unfortunately, clearly, this expression is not Lorentz invariant because of the use of the non-Lorentz invariant integration measure d3 *p*. We review below how Lorentz-invariant quantum states |*p* may be redefined so that Fourier expansions like (A29) can be made manifestly covariant.

In order to expand the quantum field into a proper continuum of plane-wave modes of the form exp <sup>±</sup>i*pμx<sup>μ</sup>*, we will need to integrate over all four-vector momenta *pμ*, i.e., perform four-dimensional integrals over **p** and *E***p** (or equivalently **k** and *ω*). However, for massive particles in general, and Proca waves in particular [94], there is a definite relation between momentum/wavevector and energy/frequency, so the integral over d3 *p* might appear at first sight to be essentially reducible to d3 *p*, while *E***p** is computed from *E*2**p** = **p**2 + *m*2. However, the problem is that the differential element d3 *p* is not Lorentz invariant, so there is a need to automatically enforce Lorentz invariance in our quantization rules. In this paper, we adopt the computationally efficient method of applying *normalizing factors* to quantum states right from the beginning. The main idea is to produce Lorentzinvariant momentum quantum states and use them to expand the quantum fields.

The essence of the method of normalizing factors is to perform a transformation of the form

$$\int\_{\mathbf{p}\in\mathbb{R}^3} \mathbf{d}^3 p \quad \underline{\qquad} \int\_{p\in\mathbb{M}^4} \mathbf{d}^4 p \tag{A30}$$

taking us from the non-Lorentz invariant space R<sup>3</sup> to Minkowski space M4, where d4 *p* is clearly invariant since *p* := *pμ* = (*p*0, **<sup>p</sup>**). However, note that the mass shell condition

$$p^2 := p^\mu p\_\mu = m^2 \tag{A31}$$

forces this integration to remain on the cone defined by equations of this type, i.e., a four-submanifold embedded into M<sup>4</sup> whose equation is (A31). Then, we may write

$$\mathbf{d}^4 p = \mathbf{d}^3 p \delta \left( p^2 - m^2 \right) \Theta \left( p^0 \right), \tag{A32}$$

where <sup>Θ</sup>(·) is the Heaviside unit step function, which is inserted by hand in order to prevent the appearance of negative energies. It can be shown that [38,76]

$$
\delta \left( p^2 - m^2 \right) \Theta \left( p^0 \right) = (1/2E\_\mathsf{P}) \delta(E - E\_\mathsf{P}) \Theta(E\_\mathsf{P}),
\tag{A33}
$$

and hence, we conclude

$$\int\_{p \in \mathbb{M}^4} \frac{\mathbf{d}^4 p}{(2\pi)^3} = \int\_{\mathbf{P} \in \mathbb{R}^3} \frac{\mathbf{d}^3 p}{(2\pi)^3} \frac{1}{2E\_\mathbf{P}} \tag{A34}$$

where the factor 1/(<sup>2</sup>*π*)<sup>3</sup> was intentionally inserted in order to make the final integral looks like an inverse Fourier transform (this inserted factor will be compensated for shortly when we define the normalizing factors of the Lorentz-invariant state |*p*). Therefore, we have managed then to reduce a four-dimensional integration in M<sup>4</sup> into a regular volume integral in R3. Note that by writing the RHS of in (A34) without *δ*(*E* − *<sup>E</sup>***p**)<sup>Θ</sup> *<sup>p</sup>*<sup>0</sup>, it is to be implicitly understood that – algorithmically speaking – every appearance of *E* or *ω* in the various possible expressions placed to the right of the integral operator sign d3 *p* should be automatically replaced by <sup>+</sup>*E***p** or +*ω***k**.

A Lorentz-invariant momentum state |*p* should be constructed from the standard three-vector momentum states |**<sup>p</sup>**, which obey the normalization rule

$$
\langle \mathbf{p} | \mathbf{p}' \rangle = \delta(\mathbf{p} - \mathbf{p}'). \tag{A35}
$$

The most obvious construction of the Lorentz-invariant momentum state then is

$$|p\rangle := (2\pi)^{3/2} (2E\_\mathbf{p})^{1/2} |\mathbf{p}\rangle. \tag{A36}$$

For example, the inner product of such states is

$$
\langle p|p'\rangle = (2\pi)^3 (2E\_\mathbf{p}) \delta(\mathbf{p} - \mathbf{p'}).\tag{A37}
$$

The definition (A36) can be readily used to yield the desired Lorentz-invariant completeness relation

$$\mathbf{1} = \int\_{\mathbf{p} \in \mathbb{R}^3} \underbrace{\mathbf{d}^3 p}\_{\text{Lorentz measure}} \frac{1}{2E\_\mathbf{P}} \mid p \rangle \langle p \vert, \tag{A.38}$$

where 1 is the unity operator. The computational utility of the method of normalizing factors stems from the fact that a completeness relation originally holding for all states {|*p*, *p* ∈ M<sup>4</sup>} in the full four-dimensional Minkowski space is now reduced to carrying out integration in the regular Euclidean space R3. However, the price to be paid for such simplification is that a momentum-dependent factor 1/2*E***p** must be inserted into the integrand of the reduced integral, which slightly complicates the evaluation of various related expressions such as probability amplitudes.

Finally, we can create the Lorentz state by applying a new normalized creation operator *b*†**p** to the ground or vacuum state |<sup>0</sup>, i.e., |*p* = *<sup>b</sup>*†**p**|<sup>0</sup> and *<sup>b</sup>***p**|<sup>0</sup> = 0. Based on (A36), we expect

$$b\_{\mathbf{p}}^{\dagger} = (2\pi)^{3/2} (2E\_{\mathbf{p}})^{1/2} a\_{\mathbf{p}\prime}^{\dagger} \qquad b\_{\mathbf{p}} = (2\pi)^{3/2} (2E\_{\mathbf{p}})^{1/2} a\_{\mathbf{p}}.\tag{A39}$$

The field expansion now, based on (A38), becomes

$$\boldsymbol{\phi}(\mathbf{x}) = \int\_{\mathbf{P} \in \mathbb{R}^3} \frac{\mathbf{d}^3 \boldsymbol{p}}{(2\pi t)^3} \frac{1}{2E\_\mathbf{P}} \left[ b\_\mathbf{P} e^{-i p\_\mu \mathbf{x}^\mu} + b\_\mathbf{P}^\dagger e^{i p\_\mu \mathbf{x}^\mu} \right],\tag{A40}$$

which after using (A39) leads to

$$\boldsymbol{\phi}(\mathbf{x}) = \int\_{\mathbf{P} \in \mathbb{R}^3} \frac{\mathbf{d}^3 \boldsymbol{p}}{(2\pi)^{3/2}} \frac{1}{(2E\_\mathbf{p})^{1/2}} \left[ a\_\mathbf{P} \boldsymbol{e}^{-i\boldsymbol{p}\_\mu \mathbf{x}^\mu} + a\_\mathbf{P}^\dagger \boldsymbol{e}^{i\boldsymbol{p}\_\mu \mathbf{x}^\mu} \right]. \tag{A41}$$

**Algorithm A1** The Canonical Quantization Algorithm (General Formulation).


### **Appendix H. On the Numerical Evaluation of the Propagator**

For the scalar massive or massless field theory considered above, the Green's Function (29) can be evaluated in closed analytical forms after regulating the integral by inserting small imaginary number i at proper locations in the integrand in order to ensure convergence (causality consideration, e.g., see [54,55]). The final expressions are given by [56]

$$G\_{\mathfrak{q}}(\mathbf{x}, \mathbf{x}') = \begin{cases} \frac{\mathrm{i}m}{8\pi} \frac{H\_{\mathfrak{i}}^{(2)}(m|\mathbf{x} - \mathbf{x}'|)}{m|\mathbf{x} - \mathbf{x}'|}, & m \neq 0, \\\frac{-1}{4\pi^2} \frac{1}{|\mathbf{x} - \mathbf{x}'|^2 - \mathrm{i}\mathbf{x}'} & m = 0. \end{cases} \tag{A42}$$

Here, the distance |*x* − *x*| is computed in Minkowski spacetime M<sup>4</sup> with the metric tensor *gμν*, i.e, the relation

$$|\mathbf{x} - \mathbf{x}'|^2 = (t - t')^2 - |\mathbf{r} - \mathbf{r}'|^2. \tag{A43}$$

The Hankel function of the second kind *H*(2) *ν* is defined as

$$H\_{\nu}^{(2)}(\mathbf{x}) := J\_{\nu}(\mathbf{x}) \pm \mathrm{i}N\_{\nu}(\mathbf{x}),\tag{A44}$$

where *Jν* and *Nν* are the Bessel and Neumann functions, respectively [95,96].

### **Appendix I. Proof of Relation**

Proof of Relation (61), when there is a source *J*(*x*), *x* ∈ *D*s, then the actual radiated (quantum) state of the q-antenna system is obtained by applying the operator

$$\int\_{\mathbf{x'} \in D\_{\mathbf{s}}} \mathbf{d}^4 \mathbf{x'} f(\mathbf{x'}) \phi^\dagger(\mathbf{x'})$$

on the ground state |0 in order to create the following one-particle excited state

$$\int\_{\mathbf{x'} \in D\_{\mathbf{s}}} \mathbf{d}^4 \mathbf{x'} f(\mathbf{x'}) \phi^\dagger(\mathbf{x'}) |0\rangle.$$

To measure the probability of having a four-momentum state |*q* in the previous excited q-antenna state, we compute the source probability amplitude *<sup>A</sup>*[*q*; *J*(*x*), *<sup>D</sup>*s] as follows:

$$\begin{split} A[q;l(\mathbf{x'}),D\_{\mathbf{s}}] \\ \vdots \\ \mathbf{f} := \langle q \vert \int\_{\mathbf{x'} \in D\_{\mathbf{s}}} \mathbf{d}^{4} \mathbf{x'} l(\mathbf{x'}) \boldsymbol{\phi}^{\dagger}(\mathbf{x'}) \vert 0 \rangle = \int\_{\mathbf{p} \in \mathbb{R}^{3}} \frac{\mathbf{d}^{3} p}{(2\pi)^{3/2} (2\omega\_{\mathbf{p}})^{1/2}} \int\_{\mathbf{x'} \in D\_{\mathbf{s}}} \mathbf{d}^{4} \mathbf{x'} l(\mathbf{x'}) \varepsilon^{\mathbf{i} p\_{\parallel} \mathbf{x'}^{\parallel p}} \langle q \vert \mathbf{p} \rangle \\ \end{split}$$
 
$$\begin{split} \mathbf{f} = \int\_{\mathbf{p} \in \mathbb{R}^{3}} \frac{\mathbf{d}^{3} p}{(2\pi)^{3/2} (2\omega\_{\mathbf{p}})^{1/2}} l(p) \delta(\mathbf{q} - \mathbf{p}) (2\pi)^{3/2} (2\omega\_{\mathbf{p}})^{1/2} = l(q) \vert\_{q^{0} = \omega\_{\mathbf{q}}} = l(\mathbf{q}, \omega\_{\mathbf{q}}) \end{split}$$

where

$$J(k) := \int\_{\mathbf{x} \in D\_{\mathbf{s}}} \mathbf{d}^{\mathbf{4}} \mathbf{x} \, J(\mathbf{x}) e^{i \mathbf{k}\_{\mu} \mathbf{x}^{\prime \mu}} = \int\_{\mathbf{x} \in \mathbb{R}^{3}} \mathbf{d}^{\mathbf{3}} \mathbf{x} \int\_{\mathbf{t} \in \mathbb{R}} \mathbf{d}t \, J(\mathbf{x}, \mathbf{t}) e^{i(\omega t - \mathbf{k} \cdot \mathbf{x})} \tag{A45}$$

is the four-dimensional (Minkowski) Fourier transform in spacetime. The second line in (A45) follows from using the momentum-space expansion of the source field *φ*(*x*) given by (A26), while in the third line, we used the normalized four-momentum state inner product formula

$$
\langle q|\mathbf{p}\rangle = (2\pi)^{3/2} (2\omega\_\mathbf{p})^{1/2} \delta(\mathbf{q} - \mathbf{p}),\tag{A46}
$$

then proceeded to evaluate the trivial resulting integral over a delta function.
