Incandescent temporal metamaterials

Regarded as a promising alternative to spatially shaping matter, time-varying media can be seized to control and manipulate wave phenomena, including thermal radiation. Here, based upon the framework of macroscopic quantum electrodynamics, we elaborate a comprehensive quantum theoretical formulation that lies the basis for investigating thermal emission effects in time-modulated media. Our theory unveils unique physical features brought about by time-varying media: nontrivial correlations between fluctuating electromagnetic currents at different frequencies and positions, thermal radiation overcoming the black-body spectrum, and quantum vacuum amplification effects at finite temperature. We illustrate how these features lead to striking phenomena and innovative thermal emitters, specifically, showing that the time-modulation releases strong field fluctuations confined within epsilon-near-zero (ENZ) bodies, and that, in turn, it enables a narrowband (partially coherent) emission spanning the whole range of wavevectors, from near to far-field regimes.

dτ ∆χ(r, t, τ )Ê(r, τ ) ·Ê(r, t), (S.I.1b) characterize, respectively, the photonic environment and the interaction term, this latter yielded by the polarization field induced by an external time-modulation [1,2]. Notice that, for clarity, we are adopting a notation showing explicitly the time dependence of the (creation and annihilation) polaritonic operators (f † and f ) [3][4][5][6]. Furthermore, it is worth noticing the mathematical character of the electric and polarization vector fields, which, within the context of macroscopic quantum electrodynamics (MQED) [6] are not functions, but quantum operators (indicated with a hat), thus bearing well known properties such as their way to be applied over quantum states, or their, in general, non-commutative character. In this regard it should be noted that the electric field operator generally reads as: As it has already been anticipated (by emphasizing the time dependence into the operators), throughout this work the analysis shall be carried out within the Heisenberg picture. It should be noted that introducing the time dependence into the states (i.e., adopting the Schrödinger picture) is not an easy task at all, since it would require for additional knowledge (or sharp assumptions) as for the dynamical regime of the system. On the contrary, under the Heisenberg picture, the dynamical behavior of the (operators characterizing the) system directly stems from the Heisenberg equation of motion: i ∂ tÔ ≡ Ô ,Ĥ . (S.II.1) Considering separately the two contributions of the Hamiltonian given in Eq. (S.I.1), it follows that: i ∂ tf (r, ω f ; t) = i ∂ tf0 (r, ω f ; t) free-evolving thus leading to: f (r, ω f ; t),Ĥ 0 =f (r, ω f ; t)Ĥ 0 −Ĥ 0f (r, ω f ; t) = ω ff (r, ω f ; t).
(S.II. 3) Proceeding similarly with the creation operator it follows that: f † (r, ω f ; t),Ĥ 0 =f † (r, ω f ; t)Ĥ 0 −Ĥ 0f † (r, ω f ; t) = − ω ff † (r, ω f ; t). Building on the above results, in this section we calculate the electric field operator making an explicit distinction between the free-evolving and the time-modulated terms. For the sake of completeness, as well as for convenience in the ensuing developments, we present both the real-valued fields in time domain, and the complex amplitudes in the frequency domain.
A. Free-evolving electric field operator The free-evolving positive-frequency electric field operator in time domain can be straightforwardly obtained just by performing the corresponding substitution of the (creation/annihilation) polaritonic operators as given in Eq. (S.II.11a) (or Eq. (S.II.11b) for the negative-frequency part), into the electric field operator given in Eq. (S.I.2a) (or Eq. (S.I.2b) for the negative-frequency part): (S.III.1) From the above expression, and noticing that we are only regarding the positive-frequency part of the field, the respective complex-like field operator in frequency domain can be directly obtained by means of the Laplace transform: where it has been used the identity, and the so-called Sokhotski-Plemelj theorem, where P stands for the Cauchy principal value, thereby accounting for the poles of causal functions. Here it is worth stressing that in order to complete the Green's function, we would have to consider the extension of the lower limit of the frequency integral to −∞. Such an approximation is commonly used (see, e.g., Refs. [1,3]), and it is often claimed that leads to some missing terms resulting from the subsequent application of the rotating wave approximation. However, in this case it should be noted that no rotating wave approximation is employed at all, so that there are no missing terms coming from such an assumption. Consequently, the second integral can be performed as follows: where the computation of the integral has been done on account of the Kramers-Kronig relations, . Therefore, putting it all together, the positive-frequency part of the free-evolving electric field operator in frequency domain reads as: B. Time-modulated electric field operator Similarly, the time-modulated electric field operator in time domain can be simply obtained by inserting the corresponding polaritonic operator into the electric field operator. Once again, focusing only into the positive-frequency part it follows that, where it has been used the completeness relation of the dyadic Green's function [3,6]: which can be derived from the Schwarz reflection principle, G * (r, r , ω) = G(r, r , −ω * ), the Lorentz reciprocity, G T (r, r , ω) = G(r , r, ω), requiring the condition that G(r, r , ω) → 0 at r → ∞ (namely, ensuring that there is no net energy transport, i.e., the time-averaged Poynting vector is to be zero for any point in space), and making use of the own definition of G(r, r , ω): In this way, the real-valued time-modulated field operator can be written as the spatial integral of a memory kernel, acting on the polarization field operator: withP(ρ; τ ) ≡ ∆χ(ρ, τ )Ê(ρ; τ ), and Just like the free-evolving contribution, the complex-like electric field operator in frequency domain can be calculated from the above expressions [Eqs. (S.III.11) and (S.III.12)] taking the Laplace transform: where it has been made the change of variablet = t − τ . Thus, the positive-frequency part of the time-modulated electric field operator in frequency domain reads as: is the Laplace transform of the polarization field operatorP(ρ; τ ) ≡ ∆χ(ρ, τ )Ê(ρ; τ ), and T . Namely, the annihilation polaritonic operator (associated to the positive-frequency component of the electric field) shall lead to a combination of the creation and the annihilation operators, being modified (actually shifted) by the time-modulation. Furthermore, in contrast to the case of a time-invariant system (e.g., that considering a quantum single emitter), there is a spatial integral that makes the kernel and the field operator tied to each other, i.e., they appear to be tightly coupled. In turn, as a consequence of the time-modulation of the system, the susceptibility function and the field operator are also coupled through the Laplace transform. At any rate, and in view of the above results, it could be inferred that, in general, for time-varying systems is impossible to express the frequency-dependent source-like (time-modulated) electric field as an algebraic product of a kernel function times a polarization field (as it actually occurs in the stationary case). This is essentially due to the coupling that exists between the medium and the external field in this kind of time-varying systems, that translates into the fact that the Laplace transform cannot be expressed in a factorized manner. Far from being a hurdle to overcome, this is indeed a distinctive feature of time-varying systems, increasing both the richness and the variety of the achievable physical effects; in particular those involving convoluted correlations that may arise from the own time-modulation of the medium. Therefore, Eq. (S.III.14), along with Eqs. (S.III.15) and (S.III.16), is arguably the most compact form in which one may express the relationship between the time-modulated electric field operator and the total field in the frequency domain. So, in order to give a more explicit result, we shall particularize the time-dependent susceptibility ∆χ to a special case, so as to be able to get an explicit expression for its associated Laplace transform.
Be that as it may, the time-modulated electric field operator,Ê T , appears both in the right-and left-hand sides of the latter expression. So, at least as a first approach, it could be useful to think on a iterative procedure in which one solves order-by-order the full expression of the time-modulated electric field operator. In this way, one could then obtain, recursively, the total electric field operator in the following way: E . According to the above result, it is clear that we can express the electric field operator as a sort of series expansion of successive higher-order terms: where: As indicated previously, to calculate the total contribution of the complex-like electric field operator in the frequency domain, we take the Laplace transform over each of the terms given above, so that: From this latter expression it can be observed that the second-order term, E (+) 2 (r; ω), gives rise to some mixing (or coupling) terms involving ancillary times appearing in the time-varying susceptibility function, ∆χ(ρ ,τ + τ ). This in turn would lead to a intertwining of Laplace transforms, hindering so the simplification of the expression to get a sort of nested (or concatenated) product of frequency-dependent functions. Still, one can take into consideration the following property of Laplace transforms relating products to convolutions: where it has been used the following definition for the Laplace transform: It is important to note that the convergence of F (ω ) should be ensured along the horizontal line Im[ω ] = ζ, that could be set to ζ → 0 + . By proceeding in this way the computation merely reduces to only one integral (eventually requiring the application of Cauchy's residue theorem), simplifying so considerably the mathematics. Then, the above results can be recast as:

IV. CURRENT DENSITY OPERATORS: DEFINITIONS AND CORRELATIONS
In this section we calculate the current density operator associated to the fields. Such a calculation can be made just by a simple comparison of the previous expressions with the following definition for the electric field operator: Similarly, for the first-order electric field, the associated current density operator reads as: Likewise, for the second-order field, and the corresponding current density operator is: These expressions allow us to perform the calculation of the corresponding nth-order correlations. To do that, we evaluate ĵ † l (ρ; ω) ·ĵ m (ρ ; ω ) th , where l + m = n, and the subscript th stands for thermal average, namely, . . . th ≡ Tr[. . .ˆ th ], withˆ th being the canonical density operator of electromagnetic fields in thermal equilibrium at temperature T (i.e., a thermal field ). Then, the zeroth-order correlation of the current density reads as: Equation (S.IV.5) is a very well-known result that provides the correlation of fluctuating currents in absence of any kind of time-modulation. Next, and from this result, one can straightforwardly calculate the first-order correlation: so that: Similarly, from the above result one could calculate the second-order correlation: In this section, we shall proceed with the calculation of the electric-field correlations of successive orders: • Zeroth-order E-field correlations: • First-order E-field correlations: where it should be noted that S 1,0 (r; ω; T ) = [S 0,1 (r; ω; T )] † . Given that the susceptibility function contains the whole information about the time-modulation, its evaluation at a constant frequency, in this case at ω = 0, gets rid of its whole dependence on the frequency and consequently, of the time-modulation. Therefore, and as a general corollary, the first-order term of the EM field correlations can be generally interpreted as an offset to the zeroth-order contribution, and, eventually, it may be totally neglected.

Angular spectrum representation of the dyadic Green's function
In order to demonstrate the usefulness of the above theoretical formulation on the fluctuation-dissipation theorem for time-varying quantum systems, here below we make use of it to tackle on a practical case in which we analyze the phenomenon of thermal emission. Specifically, we consider a flat interface separating a semi-infinite homogeneous and nonmagnetic medium (z ≤ 0), held in local thermodynamic equilibrium at a uniform temperature T , from a vacuum (z > 0). For such a system, and in general, for any translationally symmetric structure exhibiting a planar geometry (namely, slabs, interfaces, or layered media, including thin-films), one can take advantage of a powerful mathematical formalism commonly referred to as the angular spectrum representation. Such an approach allows one to express the dyadic Green's function of a system as a superposition of elementary plane waves [8,9]: Then, assuming a homogeneous, isotropic and linear * medium, it follows that [10]: Notice that, solely in the case of lossless media (i.e., those in which ε i (ω) is real), it is possible to establish the following usual correspondence: =⇒ propagating modes (far-field), which for the vacuum leads to the most standard form: Regarding this general framework, henceforth we shall only deal with two different media: • medium 1, that will be the vacuum, so that ε 1 (ω) ≡ 1, and thus having a longitudinal wavevector such that k z,1 = k 0kz,1 = k 0 κ z,1 , with κ z,1 = 1 − κ 2 R , which will be denoted by means of (x, y, z)-coordinates; • medium 2, that will be a dispersive material, so that ε 2 (ω) ≡ ε(ω), which, in general will be a complex-like function, and thus having a longitudinal wavevector k z,2 = k 0kz,2 = k 0 ε(ω)κ z,2 , with κ z,2 = 1 − κ 2 R /ε(ω) (that, due to the complex-like character of ε(ω), it will not directly fulfill the above correspondence with propagating and evanescent modes). In this case, we shall use (ρ x , ρ y , ρ z )-coordinates to denote this medium † . * Notice that, even though strictly speaking the time-modulation would break down the linear behavior of the medium, according to the Hamiltonian given in Eq. (S.I.1), we are actually assuming a linear response to relate the polarization to the external electric field. † It should be noted that, when analyzing the contributions resulting from the time-modulation we will need to deal with some auxiliary coordinates of the same medium associated to different currents; they will be denoted with the same coordinates by adding some tildes.
Likewise, and on this basis, it is worth setting down a general formulation to deal with the dyadic notation of the Green's tensors (cf. Refs. [8][9][10][11] for further details on it). Specifically, considering a plane interface system as that mentioned above, the Green's tensor relating currents in the lower half-space (medium j; with coordinates ρ z ≤ 0) to fields in the upper half-space (medium i; with coordinates z > 0), can be generally expressed as: result from the dyadic product of the corresponding polarization-vector basis, with the signs + and − being associated to modes propagating upward (↑) and downward (↓). Furthermore, is the field propagator, where it has been implicitly assumed that the medium j is underneath the medium i.
On the other hand, in the case of two points lying in the same medium, e.g., in the lower medium j, the above Green's tensor has to be recast as: are the Fresnel reflection coefficients.

Zeroth-order of thermal emission
Taking into account the previous section, the zeroth-order electric field correlation can be recast as follows: Then: where G 0,0 (r, ω) = with V being the volume of the material system (i.e., medium 2) occupying the half-space z ≤ 0.
Calculation of G0,0: We now proceed with the explicit calculation of G 0,0 . To do that, we start by substituting Eq. (S.VI.1) into (S.VI.15): where the subscript F stands for the so-called Frobenius norm, defined as: (S.VI.16) Thus: where, in accordance with the aforementioned notation, we have set the correspondenceĜ(κ x , κ y ; ω|z, ρ z ) → 2πĜ (↑/↑) z←ρz (κ R , κ ϕ ; ω|z, ρ z ), just by making a change of variables from Cartesian to cylindrical coordinates, thus leading to the factor 2π which results from the integration over the angles κ ϕ . ‡ We have also included the superscript "(↑ / ↑)" as well as the subscript "z ← ρ z " to indicate the polarization as well as the coordinates associated to the involved media. In this way, the explicit expression of the dyadic Green's function in the angular spectrum representation for this case is given by (cf. Refs. [12,13]): z←ρz (κ R , κ ϕ ; ω), the Fresnel transmission coefficients: Then, inserting Eqs. (S.VI.19a)-(S.VI.20b) into (S.VI.18), the dyadic Green's function reads as: (S.VI.22c) ‡ The change to cylindrical coordinates in the wavevector is so that κx = κ R cos κϕ, κy = κ R sin κϕ, and κz = κz, with κ 2 Thus, putting it all together it follows that: Finally: Second-order of thermal emission As shown above, the zeroth-order correlation, S 0 (r; ω; T ), leads to the ground spectra in which the time-modulation does not enter into play. In this sense, the first-order correlation S 1 (r; ω; T ) should provide for the first insights about the role of time-modulation into the system, e.g., showing up some non-local features described by means of the spatial integral. Notwithstanding, it can be observed that, at such an order, the time-modulation, characterized by means of the Laplace transform of the time-dependent susceptibility, does not account for effective frequency correlations. Essentially, this is because, in the calculation of thermal emission spectra, one is actually interested in correlations evaluated at a given frequency. For the first-order, only frequencies that are separated by the material's modulating frequency can correlate, so that, according to the result shown in Eq. (S.V.2), currents at the same frequency would not contribute effectively to the first-order spectrum. Hence, the first term in the development yielding a non-null contribution due to the time-modulation would be the second one. As pointed out previously, this term can be expressed as a sum of three contributions: S 2 (r; ω; T ) = S 2,0 (r; ω; T ) + S 0,2 (r; ω; T ) + S 1,1 (r; ω; T ), each of them displaying different behaviors depending on the specific currents involved. In this regard, correlations S 2,0 (r; ω; T ) + S 0,2 (r; ω; T ) will connect different frequencies appearing at different locations of the system, thus introducing a kind of non-local behavior that preserves the distribution (and so, the average energy) of the oscillators. On the other hand, the term S 1,1 (r; ω; T ) owns the additional feature of affecting to the photon distribution, producing a frequency-shift related to the frequency of modulation. In the following, we shall explicitly look into each of these contributions for the particular case pointed out above.
Looking at the results given in the previous section, the second-order correlations of the electric-field are given by: Likewise, the current density correlations involved in the above expressions are given by: For simplicity, as well as convenience, henceforth we shall particularize to time-harmonic modulations of the form: ∆χ(r, t) = ε 0 δχ∆χ(t) = ε 0 δχ sin (Ωt) = ε 0 δχ 2i e iΩt − e −iΩt , (S.VI.27) and whose Laplace transform is simply given by: thereby leading to: and where it has been defined ω ± ≡ ω ± Ω. Similarly, so that, , (S.VI.32) withΩ ≡ Ω/ω, and To carry on with the calculation ofḠ (ω ± ) 1,1 , first of all it is important to ensure that all of the matrices involved are properly arranged: In that way it is easy to see that the four matrices can be grouped as: , noticing that the summation over γ and δ actually refers to the customary matrix product. Inserting this latter expression into (S.VI.33a), and making use of the expansion given in Eq. (S.VI.1), it follows that: According to Eq. (S.VI.34), this latter expression has to be integrated over ρ z , ρ z , and ρ z . To do so, it should be noted the presence of the absolute value in some of the exponents so as to distinguish between the four possible cases: • Case (↑ / ↑): ρ z > ρ z and ρ z > ρ z ⇐⇒ (ρ z , ρ z ) > ρ z ; • Case (↑ / ↓): ρ z > ρ z and ρ z < ρ z ⇐⇒ ρ z > ρ z > ρ z ; • Case (↓ / ↑): ρ z < ρ z and ρ z > ρ z ⇐⇒ ρ z > ρ z > ρ z ; • Case (↓ / ↓): ρ z < ρ z and ρ z < ρ z ⇐⇒ (ρ z , ρ z ) < ρ z .
Therefore, the spatial integral in Eq. (S.VI.34) can be recast as: thereby covering the whole time-modulated medium. Obviously, all of these integrals are to be put in correspondence with the dyadic Green's functions given above. Thus, by putting it all together, the first term of Eq. (S.VI.36) can be expressed as: Proceeding similarly with the rest of contributions of Eq. (S.VI.36): Therefore, as in the previous case, Eq. (S.VI.47) can be recast as: , G (ω ± )(↓/↓) .