Broadband circularly polarized thermal radiation from magnetic Weyl semimetals

We numerically demonstrate that a planar slab made of magnetic Weyl semimetal (a class of topological materials) can emit high-purity circularly polarized (CP) thermal radiation over a broad mid- and long-wave infrared wavelength range for a significant portion of its emission solid angle. This effect fundamentally arises from the strong infrared gyrotropy or nonreciprocity of these materials which primarily depends on the momentum separation between Weyl nodes in the band structure. We clarify the dependence of this effect on the underlying physical parameters and highlight that the spectral bandwidth of CP thermal emission increases with increasing momentum separation between the Weyl nodes. We also demonstrate using recently developed thermal discrete dipole approximation (TDDA) computational method that finite-size bodies of magnetic Weyl semimetals can emit spectrally broadband CP thermal light, albeit over smaller portion of the emission solid angle compared to the planar slabs. Our work identifies unique fundamental and technological prospects of magnetic Weyl semimetals for engineering thermal radiation and designing efficient CP light sources.

In this work, we reveal that recently discovered magnetic Weyl semimetals, which are strong gyroelectric or magneto-optic materials, can emit strong circularly polarized thermal radiation which is broadband in the mid-and long-wave infrared wavelength range.
Although linearly polarized thermal radiation has been demonstrated using patterned or microstructured geometries [19,20], realizing and detecting circular polarization of thermal emission is more difficult and has been demonstrated in very few experiments using either a chiral absorber metasurface [21] or a polarization conversion metasurface [10]. Apart from these experiments and many theoretical proposals relying on nanophotonic design of reciprocal chiral absorbers [22][23][24], new approaches based on nonequilibrium coupled antennas [25] and nonreciprocal media [26][27][28] have been recently explored for CP thermal emission. However, in all of these works so far, either the purity of circular polarization is low or high-purity circular polarization is achieved over configuration-dependent narrow spectral and angular bandwidths. Here, purity can be quantified as the ratio (I (LCP) − I (RCP) )/(I (LCP) + I (RCP) ) where LCP/RCP denotes left/right circularly polarized light and I denotes the intensity. In Ref. [10], 28% purity CP thermal emission in normal direction from a metasurface is realized over 8 − 12µm range while in other works [21][22][23]25], greater than 50% purity CP thermal emission is reported but over a small wavelength band of 1 − 2µm and only along specific emission directions. In this work, we show that a rather simple geometry of a planar slab made of a representative magnetic Weyl semimetal [29], without any patterning or application of external magnetic field, can lead to greater than 50% purity CP thermal emission over a bandwidth 6 − 14µm. for around 75% of its emission solid angle.
Weyl semimetals are topological quantum materials whose nontrivial bulk band topology manifests as unique physical effects such as the topologically protected Fermi arcs, the chiral anomaly in the bulk, and large negative magnetoresistance [30,31]. The topological properties of Weyl semimetals exist only in materials which lack either time-reversal symmetry (magnetic Weyl semimetals) or inversion symmetry (noncentrosymmetric Weyl semimetals) or both. In this work, we only focus on magnetic Weyl semimetals because of their giant nonreciprocity and henceforth, we abbreviate them as MWS for brevity. Theoretically, there are many candidate MWS materials [32][33][34][35][36], and some of them such as Co 3 Sn 2 S 2 [37,38] large optical nonlinearity [41], large photovoltaic [42] and photogalvanic [43] effects, MWS are promising because of their giant nonreciprocity at infrared wavelengths which is unmatched by the existing material options of doped semiconductors in external magnetic field and conventional magnetic materials [44,45]. Because of this advantage which stems from the large intrinsic anomalous Hall effect, the nonreciprocity of MWS has been recently exploited for designing ultrasmall footprint optical isolators [46] and for optimizing thermal energy harvesting by nonreciprocal violation of directional Kirchhoff's law of thermal radiation [29,47,48]. Our work highlights yet another unique prospect of utilizing giant nonreciprocity of MWS for realizing broadband CP infrared radiation. For planar slabs, we demonstrate that the bandwidth of CP thermal light is strongly correlated with the separation of Weyl nodes in the momentum space. Figure 1 displays the planar geometry considered in this work. Using a representative Weyl semimetal analyzed in Ref. [29], the figure plots the spectra of RCP and LCP emissivities (η (RCP/LCP) ∈ [0, 1]) along the normal emission direction for momentum separations given by b = 2nm −1 (top) and b = 4nm −1 (bottom). Larger the momentum separation between the Weyl nodes, larger is the bandwidth of emitted CP thermal light. We elucidate the dependence of this result on various physical parameters describing the MWS such as Fermi energy, number of Weyl nodes, background permittivity etc., and point out the favorable parameters for achieving larger bandwidths.
We also demonstrate the circular polarization of thermal emission from finite subwavelength as well as wavelength-size bodies made of MWS. We use the thermal discrete dipole approximation (TDDA) method suitable for thermal radiation calculations from finite nonreciprocal bodies [49], and extend it to analyze the spin or circular polarization of thermal radiation from cube and rod geometries. Our computational simulations and analysis show that strong circular polarization (spin) of thermal radiation is not limited to an extended planar slab but can also be observed with finite-size objects made of MWS.
We note that chiral absorption of Weyl semimetals has been noted recently [50,51]. However, the implications for CP thermal radiation and its broadband nature have not been reported till date. A chiral absorber can emit partially CP thermal radiation based on a version of Kirchhoff's law which equates emissivity with absorptivity separately for RCP and LCP light. However, this photon-spin-resolved Kirchhoff's law is applicable only for reciprocal media [52,53]. It does not hold true for nonreciprocal media such as MWS considered here. Interestingly, nonreciprocal media follow different modified forms of spin-resolved Kirchhoff's laws [28]. The theoretical techniques and computational tools for analyzing the spin or circular polarization of thermal radiation from nonreciprocal bodies are only recently developed in the field of thermal radiation [28,54,55]. All prior inquiries concerning the topic of thermal spin photonics analyzed reciprocal media for which analytic derivations and computational simulations are simplified by invoking the principle of electromagnetic reciprocity. That is no longer the case for nonreciprocal media.
The manuscript is organized as the following. We provide the main results in section (II) followed by the conclusion (III) and the appendix (IV). In section (II), we provide the description of optical response or permittivity of MWS (II A) followed by the detailed analysis of spin or circular polarization of thermal emission from extended planar slabs (II B) and finite-size objects (II C).

A. Magnetic Weyl semimetals
There has been a significant body of research in the past decade devoted to the fundamental understanding of Weyl semimetals.
The band structure of Weyl semimetals contains an even number of nondegenerate band-touching points called as Weyl nodes which are topologically stable.
These Weyl nodes appear in pairs of two nodes of opposite chirality separated by a vector 2b in momentum space and by a scalar 2 b 0 in energy space (chiral chemical potential). The iii electromagnetic response of Weyl semimetals is described using the formalism of axion electrodynamics where these parameters appear in the form of axion term in the Lagrangian density. The axion term eventually leads to the following modified constitutive relation for the displacement field in the frequency domain [29,57]: where ε d (ω) is the frequency-dependent permittivity of the corresponding Dirac semimetal, which has b 0 = b = 0. In this work, we consider only materials where Weyl nodes have the same energy such that b 0 = 0. Throughout the work, we choose b = bk z to be along thê k z direction of the coordinate system but also mention the implications of other directionalities. Based on these considerations, the permittivity tensor of the MWS is: The presence of off-diagonal anti-symmetric entries in the permittivity tensor leads to electromagnetic nonreciprocity [58]. Such a material is called as gyrotropic medium with gyrotropy axis along z direction leading to nonzero ε xy , ε yx terms. The diagonal part of the permittivity ε d is [29,44,45]: Here ε b is the constant background permittivity, and σ is the bulk conductivity, whose dispersion is given by beyond which the band dispersion is no longer linear [44]. In the following, we first use a representative example with parameters similar to Ref. [29] as ε b = 6.2, τ = 1000 fs, ξ c = 3, b = 2 × 10 9 m −1 , ν F = 0.83 × 10 5 m/s and Fermi energy E F (T ) = 0.15 eV at T = 300 K. Then we vary these parameters within the range of reported values in the literature and analyze the circular polarization of thermal emission.

B. Extended Planar Slabs
We consider a semi-infinite half space (thick planar slab) of representative MWS (parameters noted above) at temperature T = 300K emitting thermal radiation into the vacuum half space (environment) at T 0 = 0K. The infinite transverse extent lies in the xy plane with normal to the surface along z direction. We analyze thermal emission in the direction characterized by the angles (θ, φ) where θ is the angle with z axis and φ is the angle made by the in-plane component of the emission-direction vector with x axis. The thermal radiation power emitted per unit surface area dA of the slab in the direction (θ, φ) per unit solid angle dΩ per unit wavelength dλ for a given polarization stateê is given by, 1 e hc/(λk B T ) −1 is the blackbody radiance at temperature T . It is divided by 2 in Eq.6 to account for two orthogonal polarization states separately. Since we are interested in the circular polarization of light, the polarization-dependent emissivity is calculated in the eigenbasis of right circularly polarized (RCP) and left circularly polarized (LCP) photons [28]. Based on the definitions of the polarization eigenstates in vacuum and the definition of the spin angular momentum density [59][60][61], a single RCP photon carries − angular momentum and a single LCP photon carries + angular momentum along the propagation direction (see appendix for convention details). We then define thermal photon spin denoted as S T for the extended planar slabs as: where η (LCP/RCP) denotes the emissivity of LCP/RCP light. Its dependence on (θ, φ, λ) in the above equation is omitted for brevity. This quantity is similar to the circular dichroism parameter considered in circular dichroism spectroscopy [62,63] or the well-known Stokes S 3 parameter [64]. It is also same as the purity of circular polarization described in the introduction. S T furthermore quantifies the net longitudinal spin angular momentum per emitted thermal photon. It can be measured experimentally using quarter wave plates and detectors oriented in the given emission direction as done in previous experiments [10]. The spin-resolved emissivities η (LCP/RCP) in the above expressions for nonreciprocal planar media have been derived previously in Ref. [28], and are separately discussed in the appendix (Section IV) in the context of this work. In direction (θ = 0, φ = 0) as a function of wavelength. A large separation between these spectra leads to |S T | ∼ 1 (black line) over a broad range of wavelengths. We note that |S T | is reasonably large beyond λ ∼ 16µm as evident from the top figure. However, total emissivity given by η T = (η (LCP) + η (RCP) )/2 is extremely small at these higher wavelengths. Since the resulting thermal radiation power is extremely small, it will be challenging to detect CP light in potential experiments. Therefore, throughout this work, we only focus on the range of wavelengths where both S T and η T are reasonably large. The gyrotropy or the strength of nonreciprocity as quantified by the ratio |ε xy |/|ε xx | 1 over the same wavelength range is displayed by the middle figure. Both top and middle figures reveal the connection between strong gyrotropy and strong CP thermal radiation. The bottom figure plots S T for all emission angles θ ∈ [0, π/2] for the same wavelength range. Because of the rotational symmetry due to perpendicular gyrotropy axis, these plots are the same for all in-plane angles φ ∈ [0, 2π]. The contour plot demonstrates that high purity (|S T | ∼ 1) CP light is emitted over not only a broad spectral bandwidth but also a broad range of emission angles (θ).
Based on the contour plot and the top figure, broadband high-purity CP thermal emission characterized by |S T | > 0.5 and η T > 0.5, is realized over almost 75% of the emission solid angle of 2π in the wavelength range of 6µm to 14µm for this representative MWS planar slab.
In Fig. 2b, we plot the same features for another well-studied nonreciprocal system of a doped InSb slab in magnetic field. For B = 1T magnetic field applied perpendicular to InSb slab of large doping concentration n = 10 17 cm −3 (permittivity response taken from previous works [55]), we find that the emissivities of RCP and LCP light differ noticeably as shown in the top figure. Again, the connection with strong gyrotropy over the same wavelength range is evident from the middle figure. The broadband nature of S T (both emission directions and wavelengths) is depicted by the contour plot in the bottom figure. However, we note that this high purity thermal spin (|S T | ∼ 1) for InSb slab is obtained at wavelengths in the range of ≈ 40 − 100µm. In typical experiments or practical applications, the operating temperatures can be considered to be in the range of T = 300 to 1000K. For such temperatures, the blackbody distribution v (I bb (λ, T )) peaks at wavelengths in the range of 3 − 12µm and decreases exponentially at larger wavelengths. The resulting low thermal emission power in addition to requirement of strong external magnetic field and limited availability of low-loss quarter wave plates over relevant wavelengths (needed for detection of nonzero S T ) are important requirements that make such doped semiconductors rather inconvenient for implementing the functionality of broadband CP thermal emission. These considerations are quite important from the perspective of high-precision experimental measurements of thermal radiation [65]. Other nonreciprocal media like ferroand ferrimagnets exhibit large magnitude of ε xy in the wavelength range of 3 − 20µm relevant for thermal radiation engineering. However, the gyrotropy remains small (|ε xy |/|ε xx | 1) because of large negative permittivity (Re{ε xx } 0). Furthermore, these materials also exhibit high reflectivity and low surface emissivity causing low thermal emission power. One can consider enhancing the gyrotropy in the relevant wavelength range based on patterning or resonant systems. However, such resonance-based approaches will necessarily lead to CP emission over narrow spectral and angular bandwidths. Optimization or inverse design of quarter wave metasurfaces (linear to CP conversion) or reciprocal chiral absorbers may be explored to obtain large bandwidths but no such work has been done so far. To the best of our knowledge, the bandwidths reported here using a representative MWS have not been realized before either theoretically or experimentally. Therefore, we believe that the planar slabs of MWS are uniquely capable of enabling the functionality of broadband CP thermal infrared radiation without patterning or application of external magnetic field.
We now tune the underlying parameters in the model described in section II A and analyze the behavior of thermal photon spin S T to point out favorable parameters for realizing broadband CP infrared radiation. We first vary the Weyl nodes separation b in the above model alongk z from 1 nm −1 to 6 nm −1 . Figure 3(a) shows the spectrum of spin S T (solid lines) as well as of total emissivity η T (dashed lines) for increasing values of b. The same convention for S T and η T is then employed for all figures in Fig. 3. Since the strength of nonreciprocity is directly proportional to Weyl nodes separation, the broadband nature of S T is further improved as the separation b is increased. For instance, both |S T | ∼ 0.5 and reasonable emissivity η T 0.25 are achieved over 4µm to 28µm wavelength range when b = 6nm −1 (violet lines). Figure 3(b) shows the dependence on Fermi energy by considering the values of E F = 0.15, 0.1, 0.07eV some of which are reported previously in [45,66,67]. For these plots, the Weyl nodes separation is b = 2nm −1 . As evident, with reduced Fermi energy, the bandwidth of CP thermal radiation increases and the spectrum shifts towards larger wavelengths because of smaller photon energies associated with the bands near low Fermi energy.
In the model of a representative MWS considered above, we used a simplified system containing only 2 Weyl nodes in the electronic band structure. This parameter, denoted as g, influences the diagonal permittivity of the Weyl semimetal through Eq. 5. Figure 3(c) demonstrates the impact of g on the thermal photon spin S T . The blue, red, yellow plots corresponding to g = 2, 12, 24 respectively are obtained when b = 2nm −1 and E F = 0.15eV. The resulting trend indicates that the bandwidth is reduced as the number of Weyl nodes is increased. The values of g = 12, 24 have been previously noted in Ref. [33,67]. In order to increase the bandwidth for a given number of Weyl nodes, we can either increase b or decrease E F as illustrated in Figs. 3(a,b). For example, comparison of violet (g = 24, E F = 0.1eV) and yellow (g = 24, E F = 0.15eV) lines in Fig.3(c) indicates the feasibility of broadband CP light with larger number of Weyl nodes by reducing E F . Figure 3(d) shows the dependence of the spectrum on the background permittivity denoted as ε b . As evident, the spectrum slightly shifts towards larger wavelengths for increasing values of ε b . Similar frequency shifts are observed for variations of reasonable magnitude in the other parameters in the model. For all results above, the temperature of the planar slab is assumed to be held at T = 300K. When the temperature is increased to T = 400K, small frequency shifts in the spectra are observed similar to those noted for total emissivity spectra in Ref. [29]. Based on the numerical simulations, we find that the parameters b, E F , g are more impactful than other parameters (ε b , ν F , τ, T ) considered in the permittivity model of MWS.
We also note that S T is quite small for the Voigt geometry in which the gyrotropy axis of MWS is parallel to the surface of the planar slab. This result was noted for InSb slab in Ref. [28] and was verified numerically for Weyl semimetal in this work. While Voigt geometry (parallel gyrotropy) is required for violation of directional Kirchhoff's law in other recent works [29,48], Faraday geometry (perpendicular gyrotropy) is found to be more suitable for realizing strong CP thermal emission. Indeed, the Voigt geometry leads to very small thermal spin for the same magnitude of gyrotropy of the MWS in the Faraday geometry. Another important aspect concerns the partial linear polarization of thermal emission because of gyrotropy of the Weyl semimetal. Similar to the analysis of emissivities in the eigenbasis of RCP and LCP light, we can derive the emissivities in the eigenbasis of transverse electric (s) and transverse magnetic (p) polarized plane waves or in the eigenbasis of s + p and s − p polarized plane waves. We note that η (s) = η (p) and η (s+p) = η (s−p) for the normal emission direction (θ = 0) because of the rotational symmetry of the Faraday geometry considered above.  angles, we find that the differences in these emissivities (|η (s) − η (p) |, |η (s+p) − η (s−p) |) are very small compared to |η (LCP) − η (RCP) |. Thus, the gyrotropy or nonreciprocity of the Weyl semimetal in Faraday geometry is primarily connected with the spin angular momentum or circular polarization of emitted thermal photons. Because of these reasons, the analysis of nontrivial partial linear polarization for θ = 0 emission directions in Voigt and Faraday geometries, which is rather straightforward, is not provided in this work.

C. Finite subwavelength and wavelength-scale objects
In addition to extended planar slabs, we now show that the finite objects made of a magnetic Weyl semimetal can also emit broadband CP thermal radiation in specific emission directions. Again, there is no local formulation of Kirchhoff's law [68] which can simplify the calculation of emission by invoking electromagnetic reciprocity. Our system is nonreciprocal and the numerical calculation of circular polarization of thermal emission from these nonreciprocal bodies is more involved than that of reciprocal bodies. One numerical approach of analyzing thermal radiation from finite-size bodies is thermal discrete dipole approximation abbreviated as TDDA [49,69,70]. Within TDDA, the finite body is subdivided into tiny subvolumes comprising of thermally fluctuating currents whose temperature-dependent strength is governed by the fluctuation-dissipation theorem (FDT). This approach was utilized to analyze radiative heat transfer in nonreciprocal systems [49,54], and was recently extended to analyze radiative angular momentum transfer [71]. Using these recent developments, we calculate the thermal photon spin S T for thermal vii emission from a finite-size body as: where W (ω) denotes the spectral energy density and the vector S(ω) denotes the spectral spin angular momentum density [59][60][61] at a given frequency ω.r denotes the emission direction. This definition of spin in the context of thermal radiation has been previously explained in Ref. [25] and it quantifies the net spin angular momentum per emitted photon in the given emission direction (r). In other words, the same quantity from the previous section is analyzed here. The quantities E(r, ω) ⊗ E * (r, ω) and H(r, ω) ⊗ H * (r, ω) are computed using the aformentioned TDDA approach applicable for finite-size nonreciprocal objects. From the perspective of potential experiments, it is also important to calculate the net emission power which has been described in several works before [49]. In the following, we plot the spectra of spin S T as well as the thermal emission power.
We first consider a cube made of MWS of length 1µm along each side. The temperature of the cube is T = 300K and that of its surrounding environment or vacuum is T 0 = 0K. We use the same representative MWS as used for the planar slab considered in Fig. 2 with gyrotropy axis alongẑ direction. Using TDDA approach by discretizing the cube into 12 × 12 × 12 subvolumes, we calculate S T in different emission directions in the far-field (over a sphere of radius 10 5 µm.). Numerical convergence in terms of less than 1% relative change is observed upon changing the discretization scheme from 10 × 10 × 10 subvolumes to 12 × 12 × 12 subvolumes. As shown in the left figure of Fig. 4(a), S T is plotted for emission directions characterized by θ ∈ [0, π]. The dependence on the angle φ in the xy-plane is negligible because of the subwavelength size of the cube and the gyrotropy axis being perpendicular to the xy-plane. The spectrum of S T reveals the broadband partially CP thermal emission with S T being larger for smaller values of θ. At θ = π/2, S T = 0 at all wavelengths and for θ > π/2, all the plots are flipped in sign (S T (π − θ) = −S T (θ)). The overall thermal emission power spectral density (integrated over all emission directions) is plotted in the right figure showing the wavelength range of interest from the perspective of potential experiments. In the wavelength range of interest between 5µm to 15µm where thermal emission power is dominant, |S T | is close to one for a broad range of wavelengths in the emission direction characterized by θ close to either 0 or π values (parallel and antiparallel gyrotropy directions). We also consider a rod of dimensions 1µm ×1µm ×10µm with its longest side oriented alongẑ direction. Numerical convergence of TDDA calculations is ensured in terms of less than 1% relative change upon changing the discretization scheme from 8 × 8 × 80 subvolumes to 10×10×100 subvolumes. The resulting spectrum of spin S T calculated in the far field is plotted in the left figure of Fig. 4(b) and the total thermal emission power spectral density is plotted in the right figure. The results are qualitatively similar to those displayed in Fig. 4(a). Both figures prove the feasibility of broadband CP radiation from finite-size bodies of Weyl semimetals albeit with a smaller range of favorable emission directions compared to the extended planar slabs.

III. CONCLUSION
We numerically demonstrated that planar slabs of magnetic Weyl semimetals (MWS) can emit high-purity CP thermal radiation (S T ∼ 1) over a broad wavelength range in mid-IR and long-wave-IR spectral region, and for a significant portion of the emission solid angle. For a representative MWS with parameters similar to Ref. [29], strong circular polarization or thermal photon spin (S T > 0.5) along with reasonably large emissivity (η T > 0.3) is realized over 6µm to 14µm wavelength range and 75% portion of the emission solid angle of 2π for the planar slab. Based on theoretically well-established model of MWS, we analyzed the dependence of this result on underlying physical parameters and showed that even larger bandwidths can be achieved by tuning the favorable parameters.
Specifically, larger separation between Weyl nodes (equivalent to stronger gyrotropy), smaller Fermi energy, and smaller number of Weyl nodes are favorable amongst other parameters for achieving larger bandwidths. We further demonstrated using TDDA computations that finite-size bodies made of MWS also lead to spectrally broadband CP thermal radiation albeit with geometry-dependent direction-specificity or small angular bandwidths compared to extended planar slabs.
We note that no such functionality of highly efficient broadband CP radiation at mid-IR and long-wave-IR wavelengths currently exists. It is difficult to realize simultaneously large spectral and angular bandwidths of high purity CP infrared emission reported in this work by employing other existing approaches. In other words, MWS appears to exhibit a unique capability to solve this technological design problem. There has been an extensive work in the past decade on Weyl semimetals with many candidates currently being investigated in experiments.

FUNDING
We acknowledge funding from DARPA Nascent Light-Matter Interactions program.

DISCLOSURES
The authors declare no conflicts of interest.

IV. APPENDIX
We provide the summary of theoretical derivations and computational methods here.

A. Spin resolved emissivity for extended planar slabs
The polarization-resolved thermal radiation power from an extended planar slab is given by Eq. 6 in the main text. This can be derived within the radiometry framework using detailed balance of energy and momenta or within the scattering formulation of fluctuational electrodynamics. The derivation for nonreciprocal media is discussed in detail in Ref. [28]. Here, we summarize the important details related to conventions of RCP and LCP light.
In the planar geometry, the emission direction is characterized by the angles θ, φ made by the propagation wavevectork.
It is straightforward to obtain the associated transverse electric (s) and transverse magnetic (p) polarization unit vectors for the plane of incidence ix spanned byk andẑ (normal to the slab surface). These eigenvectors are given below: The RCP eigenvector isê s + iê p (also denoted as +) and the LCP eigenvector isê s − iê p (also denoted as −). From the perspective of an observer/detector in the far-field, the incoming RCP photons have a clockwise rotation of E-field and the incoming LCP photons have an anti-clockwise rotation of E-field. The spin angular momentum per photon along the propagation direction is for LCP photons and − for RCP photons [64]. The emissivities in the eigenbasis of RCP and LCP photons are: where R (ij) (ω, θ, φ) for i, j ∈ {+, −} denotes the polarization interconversion reflectance for light of angular frequency ω incident in the direction characterized by the angles (θ, φ). These reflectances depend on the associated Fresnel reflection coefficients inê s ,ê p basis and given by the following expressions: R (−+/+−) = |(r ss − r pp ) ± i (r sp + r ps )| 2 /4.
where we have omitted (ω, θ, φ)-dependence for brevity since it is the same on both sides. r jk (ω, θ, φ) denotes the amplitude of j-polarized reflected light due to unit-amplitude incident k-polarized light of frequency ω.
These reflection coefficients are computed by solving the boundary conditions. A computational tool to obtain these coefficients for generic bianisotropic media is available on GitHub https://github.com/chinmayCK/Fresnel under the MIT license.

B. Thermal discrete dipole approximation approach
We employ thermal discrete dipole approximation method recently used for analyzing thermal emission from nonreciprocal finite bodies [49]. Here we clarify the assumptions and details in TDDA formalism specific to this work. Within TDDA formalism employed in this work, we discretize the finite-size objects into cubic sub-volumes containing point dipoles at their respective centers. Throughout this work, we assume the environment to be at 0K. The electric field outside the object generated by discrete fluctuating point dipoles is expressed using dyadic Green's tensors as: whereP is a column vector of length 3N consisting of dipoles p 1 , . . . , p N , and G EE = Ĝ EE (r, r 1 ), · · · ,Ĝ EE (r, r N ) is the row vector of (17) where R = r − r , R = |r − r |, and ⊗ denotes the outer product. Here the bar notation ofP and other following variables indicates that the quantities are located or evaluated inside the object. The discrete dipoles inside the object denoted as p i consist of fluctuating part p  These equations in the matrix form are as the following: P =P (fl) +P (ind) , This forms the discrete dipole approximation (DDA) equation which is to be solved to obtainP required to calculate thermal emission quantities. The fluctuating point dipoles p (fl) i appearing in the above equations satisfy the fluctuation-dissipation theorem (FDT). For our calculations, the FDT leads to following expression for correlations in the matrix form: P (fl) (ω)P (fl) † (ω ) = 2π ε 0 δ(ω − ω )(1 + 2n B (ω, T ))χ (19) where n B (ω, T ) = 1/(e ω/k B T − 1) is the Bose-Einstein distribution function, andχ = diag(χ 1 , . . . ,χ N ),χ i = x χ is a tensor combining the imaginary part of the polarizability tensor and the radiative correction, Andα = ( 1 Vp (L p + [ε p − I] −1 ) − i k 3 0 6π I) −1 is the polarizability tensor for cubic subvolume containng a material characterized by a generic permittivity tensor ε p which is not necessarily diagonal as in the case of isotropic materials. V p is the volume of the cubic element. k 0 = ω/c where c is speed of light andL p = (1/3)I.
By solving this DDA equation given by Eq.18, we get P =T −1P(fl) whereT = I − k 2 0ᾱ dḠ EE . Therefore, from which we compute the correlation E(r, ω) ⊗ E * (r, ω) , which is Similarly, the magnetic field correlations H(r, ω) ⊗ H * (r, ω) are obtained by replacing the electric-electric dyadic Green's tensorĜ EE by magnetic-electric tensor: Similarly, the correlations E(r, ω)⊗H * (r, ω) required for calculation of Poynting flux or radiative heat transfer can be obtained within the TDDA formalism, which has been discussed in several works before [49,69,70,72]. We employ the approach taken in Ref. [49] in this work for calculating the angle-integrated thermal emission spectral density in the far field, plotted in Fig.4