Thermal spin photonics in the near-field of nonreciprocal media

The interplay of spin angular momentum and thermal radiation is a frontier area of interest to nanophotonics as well as topological physics. Here, we show that a thick planar slab of a nonreciprocal material, despite being at thermal equilibrium with its environment, can exhibit nonzero photon spin angular momentum and nonzero radiative heat flux in its vicinity. We identify them as the persistent thermal photon spin (PTPS) and the persistent planar heat current (PPHC) respectively. With a practical example system, we reveal that the fundamental origin of these phenomena is connected to spin-momentum locking of thermally excited evanescent waves. We also discover spin magnetic moment of surface polaritons in nonreciprocal photonics that further clarifies these features. We then propose a novel thermal photonic imaging experiment based on Brownian motion that allows one to witness these surprising features by directly looking at them using a lab microscope. We further demonstrate the universal behavior of these near-field thermal radiation phenomena through a comprehensive analysis of gyroelectric, gyromagnetic and magneto-electric nonreciprocal materials. Together, these results expose a surprisingly little explored research area of thermal spin photonics with prospects for new avenues related to non-Hermitian topological photonics and radiative heat transport.


I. INTRODUCTION
Thermal spin photonics merges the fields of the thermal radiation and the spin angular momentum of light. Thermal radiation plays an important role in energy-conversion and renewable technologies [1,2] while the spin angular momentum property of light is fundamentally relevant in the context of spin-controlled nanophotonics [3,4], chiral quantum optics [5] and spintronics [6]. Despite extensive work in the past few decades, there has been very little overlap between these two areas.
Important developments in this field include spin-polarized (circularly polarized) far-field thermal radiation from chiral absorbers [7][8][9] and the definition of the degree of polarization in the thermal near-field [10] of reciprocal media. In stark contrast, the primary aim of this work is to explore thermal spin photonics (spin-related thermal radiation phenomena) in the near-field of non-reciprocal media.
Our work utilizes fluctuational electrodynamics and is fundamentally beyond the regime of Kirchhoff's laws which is valid only for far-field thermal emission from bodies at equilibrium.
One striking example where spin angular momentum of thermal radiation is not captured by Kirchhoff's laws, is circularly polarized thermal emission from coupled non-equilibrium antennas demonstrated in our recent work [11]. This approach of exploiting interacting non-equilibrium bodies is fundamentally unrelated to conventional approaches of achieving spin angular momentum of light based on either polarization conversion or structural chirality. Practically, this non-equilibrium mechanism enables temperature-based reconfigurability of the spin state of emitted thermal radiation. Our current work deals with bodies at thermal equilibrium with their surroundings, and reveals surprising spin angular momentum features in their near-field arising in presence of nonreciprocity.
To show the universal nature of these non-reciprocal thermal spin photonic effects, we develop a framework to analyze equilibrium thermal-radiation properties of a planar slab of a generic bianisotropic material described by arbitrary permittivity (ε), permeability (µ) and magneto-electric susceptibilities (ξ, ζ). Such a material is reciprocal if the material properties satisfy, It is nonreciprocal if any one of these conditions is violated. With fluctuational electrodynamic analysis, we show that a nonreciprocal planar slab at thermal equilibrium with its environment can exhibit nonzero spin angular momentum of thermal radiation in its near-field. We identify it as the persistent thermal photon spin (PTPS) because it exists without any temperature difference analogous to well-known persistent electronic charge current [12,13] that exists without any voltage difference. The PTPS is also accompanied by locally nonzero radiative heat flux parallel to the surface which we call as the persistent planar heat current (PPHC).
We reveal that the spin-momentum locking [14][15][16] of thermally excited evanescent waves, plays a fundamental role in facilitating both these phenomena with a practical example system.
Our work thus provides the first generalization of the spin-momentum locking of light well-known in topological photonics and atomic physics to thermally excited waves. We consider a doped Indium Antimonide (InSb) slab at room temperature with an arbitrarily directed magnetic field. The thermally excited surface plasmon polariton supported by InSb slab has transverse spin locked to its momentum [14]. Our calculations reveal that ii the spin-momentum locked polariton (electromagnetic wave) also carries spin magnetic moment which leads to polaritonic energy/frequency shift through Zeeman type interaction with the applied magnetic field. For InSb sample with doping concentration ∼ 10 17 cm −3 , the polaritonic spin magnetic moment is found to be around 10µ B where µ B is Bohr magneton. The polaritonic magnetic moment depends asymmetrically on the momentum for forward and backward propagating polaritons leading to asymmetric energy shifts. This clarifies the fundamental origin of PTPS and PPHC, resulting from asymmetric contributions of forward and backward propagating evanescent waves.
Detecting thermal radiation effects of non-reciprocal media is an open challenge. We note however, that our discovered effects PTPS and PPHC are significantly enhanced in the near-field due to a large density of thermally excited evanescent states. In particular, we show the striking result that at a distance d 0.5µm from the slab surface, the magnitude of PTPS exceeds the spin angular momentum density contained in the laser light carrying typical power of ∼ 1mW. This immediately motivates experimental validation of our predicted effects by probing optical forces and torques on small absorptive particles in the thermal near-field of an InSb slab. We predict that the Brownian movement of these particles will be sufficiently influenced by the additional thermal spin photonic forces and it can be directly viewed using a lab microscope.
We further demonstrate the universal behavior of both these near-field thermal radiation phenomena with a comprehensive analysis of the key classes of nonreciprocal media namely, gyroelectric (ε = ε T ), gyromagnetic (µ = µ T ) and magneto-electric (ξ = −ζ T ) materials. The general analysis describes the origin and the nature of these features for any given material type and further reveals that a nonreciprocal material is necessary but not sufficient to observe PTPS and PPHC.
Our work advances science in multiple directions. It makes a new fundamental connection between the spin-momentum locking of evanescent waves [14][15][16] and the radiative heat transfer. The spin-momentum locking of thermally excited waves opens a new degree of freedom for directional heat transport at the nanoscale. The spin magnetic moment of gyrotropic surface polaritons (which contain both s-and p-polarized waves) invites related studies of spin-dependent quantum plasmonics and spin-quantization. We also address the experimental detection of the persistent thermal photonic phenomena which is not addressed by previous works [17][18][19]. It is important for thermodynamic revalidation of fundamental understanding of nonreciprocal systems and also because there is no experiment till date probing the intriguing effects of nonreciprocity on thermal radiation. While thermal photonics is so far limited to isotropic, anisotropic and gyroelectric materials [20,21], we provide a universal description for all material types to motivate similar studies of thermal-radiation phenomena with largely unexplored material types such as topological insulators, multiferroic and mangeto-electric materials. We note that the theoretical framework and the tools employed here will be useful for studying not only fluctuational [20][21][22][23][24][25] but also quantum [26][27][28][29][30] electrodynamic effects by using generic, bianisotropic materials.

II. RESULTS
Theory. We consider a planar geometry shown in figure 1 comprising of a semi-infinite half-space of generic homogeneous material, interfacing with semi-infinite vacuum half-space at z = 0.
We focus on the thermal radiation on the vacuum side of this geometry where the physical quantities such as energy density W , Poynting flux P and spin angular momentum density S are well-defined [31,32] and measurable in suitable experiments [33,34]: where r denotes the position vector and ... denotes the thermodynamic ensemble average. The spin angular momentum density (3) has so far been studied primarily for non-thermal light [31,32], where it leads to proportionate optical torque on small, absorptive particles [35,36]. We have generalized it here and in our recent work [11] to thermally generated electromagnetic fields in vacuum.
We calculate both electric and magnetic type thermal spin angular momentum density given by S (E) and S (H) respectively. Throughout the manuscript, all quantities are described in SI units and the dependence on frequency ω (such as E(ω), H(ω)) is suppressed assuming e −iωt time dependence in Maxwell's equations. The above quantities Q = {W, P, S} are to be integrated over frequency to obtain the total densities/flux rates asQ = ∞ −∞ dω 2π Q(ω)dω. Keeping in mind the future explorations using generic, bianisotropic materials, we prefer to use vector potential in Landau gauge to obtain the electromagnetic fields (E = iωA, H = ∇ × A). The electromagnetic field correlations required for calculation of densities and flux rates above are obtained from the vector potential correlations. These correlations evaluated at two spatial points r 1 , r 2 are expressed in the matrix form as A(r 1 ) ⊗ A * (r 2 ) = Geometry. We analyze near-field thermal radiation properties namely energy density W (r), Poynting flux P(r) and photon spin angular momentum density S(r) for a generic bianisotropic planar slab characterized by permittivity ε, permeability µ and magneto-electric coupling tensors ξ, ζ. The yellow ovals indicate underlying fluctuating dipoles that emit thermal radiation. We show that thermal photon spin density S(r) and heat flux P(r) can be nonzero for a nonreciprocal material despite thermal equilibrium between vacuum and planar slab.
Here Θ(ω, T ) = ω/2 + ω/[exp( ω/k B T ) − 1] is the average thermal energy of the harmonic oscillator of frequency ω at temperature T . The Green's tensor G(r 1 , r 2 ) relates the vector potential A(r 1 ) to all the source currents J(r 2 ) such that A(r 1 ) = We derive the Green's function given in the methods section for a planar slab of a generic, bianisotropic medium (see supplement for derivation).
Finally, using the Green's function and vector potential correlations above, we obtain the electromagnetic field correlations at two spatial points r 1 and r 2 in vacuum: where ∇ rj × for j = [1, 2] is the differential curl operator.
The densities and flux rates are then calculated using definitions (1),(2) and (3) using above correlations with r 1 = r 2 = r and making use of the fluctuation-dissipation relation given by Eq.(4). Possibility of observing nonzero spin angular momentum and heat flux despite thermal equilibrium.
We now show that a nonreciprocal medium can lead to nonzero spin angular momentum density (3) and nonzero Poynting flux (2) in its thermal near-field. We make use of insightful expressions derived here and time-reversal symmetry arguments, but we do not refer to any specific material in this section.
The electromagnetic waves in this planar geometry are characterized by their in-plane conserved propagation wavevector k. Poynting flux and spin angular momentum of each wave are P(k) and S(k) respectively. It follows from time-reversal symmetry that heat and angular momentum associated with thermally excited k-waves are negated by heat and angular momentum carried by −k-waves (P(k) = −P(−k), S(k) = −S(−k)), resulting into zero flux rates at thermal equilibrium. Because of violation of the time-reversal symmetry for nonreciprocal media, this is no longer true and one can expect to see nonzero flux rates in the absence of cancellation. We identify the resulting nonzero spin angular momentum density as the persistent thermal photon spin (PTPS) and nonzero heat flux as the persistent planar heat current (PPHC). Although this analysis proves that nonreciprocity is necessary to observe PTPS and PPHC, full fluctuational electrodynamic calculations below, confirm and reveal much more, including the result that nonreciprocity is not a sufficient condition.
Before demonstrating the full calculations, we can make some general comments regarding heat flux and thermal photon spin perpendicular to the slab for which semi-analytic expressions are insightful (see supplement for their derivation). The near-field Poynting flux inê z direction, for any material at thermal equilibrium with vacuum half-space. Similarly, the electric and magnetic parts of spin angular momentum density alongê z direction are: It follows that the total perpendicular thermal spin , is zero for any material at thermal equilibrium.
Here r ss , r pp , r sp , r pp are the Fresnel reflection coefficients for light incident on the planar iv slab having perpendicular wavevector k z , (conserved) parallel wavevector k , and making an azimuthal angle φ with x-axis of the geometry. k 0 = ω/c = k 2 z + k 2 is the vacuum wavevector (see methods and supplement). For reciprocal media, (r sp + r ps ) = 0 [37] and the electric-and magnetic-type persistent thermal photon spin (PTPS) are separately zero. On the other hand, this condition is not necessarily true for nonreciprocal media and interestingly, even though electric-and magnetic-type PTPS can be separately nonzero, total PTPS perpendicular to slab is always zero. We note that the semi-analytic expressions for PTPS and PPHC parallel to the slab and given in the supplement are not conducive for such general insights but full calculations of examples below reveal their existence and nature.
It is important to point out here that the intuition confounding presence of non-zero heat current at thermal equilibrium does not lead to thermodynamic contradictions. In particular, because of the nonzero heat flow parallel to the surface, it could be expected that one end will be hotter than the other end. However, given the infinite transverse extent of the system considered above, there is no end that can be heated or cooled [38]. On the other hand, since the two distinct half-spaces are separated by a well-defined interface, no macroscopic flux rates can exist across the boundary by definition of thermal equilibrium between the half-spaces. The fluctuational electrodynamic theory produces a consistent result above that the Poynting flux (P z ) and the total spin angular momentum density (S z ) perpendicular to the surface are identically zero for any material. For finite-size nonreciprocal systems such as finite planar slabs or other geometries (cylinders, cubes) having well-defined edges, it follows from the energy conservation under global thermal equilibrium that, the energy exchange of any finite subvolume V of the system with the rest of the system is zero i.e. ∂V P.dA = 0 where ∂V denotes the surface and dA is the differential area vector. It then follows from the divergence theorem that the Poynting flux P is divergence-free everywhere (∇ · P = 0). This means that there are no sources or sinks for Poynting vector lines and they form closed loops. Therefore, the persistent current in the near-field of a finite-size nonreciprocal system will flow around the edges and form a closed loop, conserving energy globally. While rigorous demonstration of fluctuational electrodynamic confirmation in arbitrary geometries is challenging, this is evident for finite spherical systems analyzed in ref. [17,18]. We also note that our analysis based on the assumption of infinite transverse extent is suitable for describing a real (finite) nonreciprocal planar slab if thermal equilibrium is achieved over transverse dimensions much larger than the wavelengths associated with persistent features. In that case, the local persistent features can be analyzed within the present theory by ignoring the edge effects. Since this situation is quite realistic as we also describe in our experimental proposal further below, our fluctuational electrodynamic analysis is adequate. Also, we remark that the planar geometry of infinite transverse extent is a reasonable theoretical approximation of planar slabs. It is well-known and extensively used in the context of closely related topics of Casimir force [39] and near-field radiative heat transport [40] between planar bodies. In the following, we consider a practical example system and demonstrate these effects and further clarify their fundamental connection with the spin-momentum locking of evanescent waves.
Practical Example of InSb slab. We consider doped Indium Antimonide (InSb) slab which has been most widely studied in context of coupled magneto-plasmon surface polaritons [41,42] and whose material permittivity dispersion has been well-characterized experimentally [43][44][45]. For the sake of completeness of our study, we extend the known permittivity model in ref. [44,45] to the case of an arbitrarily oriented magnetic field in our geometry and obtain the semi-analytic form of permittivity given This is obtained from an extended Lorentz oscillator model in which bound/free electrons (charge q, effective mass m f , position r) are described as mechanical oscillators, that further experience additional Lorentz force (qṙ × B) in presence of applied magnetic field B. Each ω cj for j = [x, y, z] describes the cyclotron frequency inê j direction given by ω cj = q(B ·ê j )/m f . In this model, ε ∞ is the high-frequency dielectric constant, ω L is the longitudinal optical phonon frequency, ω T is the transverse optical phonon frequency, ω p = nq 2 m f ε∞ε0 is the plasma frequency of free carriers of density n and effective mass m f where q is electron charge and ε 0 is vacuum permittivity. Γ denotes the optical phonon damping constant while γ is the free-carrier damping constant. All parameters are obtained from Ref. [45] for InSb sample of doping density n = 10 17 cm −3 . ε ∞ = 15.7, ω L = 3.62 × 10 13 rad/s, ω T = 3.39 × 10 13 rad/s, ω p = 3.14 × 10 13 rad/s, Γ = 5.65 × 10 11 rad/s, γ = 3.39 × 10 12 rad/s, m f = 0.022m e where m e = 9.1094 × 10 −31 kg is electron mass. Because of the anti-symmetric (gyroelectric-type) permittivity tensor (ε = −ε T ) of InSb in presence of magnetic field, it is nonreciprocal and can lead to persistent features in its thermal near-field. For the above paramenters, InSb slab supports surface plasmon polaritons (SPPs) at ω ∼ 3.9 × 10 13 rad/s v ф = π/2 ф = 0 ф = -π/2 ω ( ×10 13 rad/s) 3 (∼ 48µm) localized close to the interface with vacuum. Because of their significant contribution to the near-field thermal radiation, we analyze these polaritons and also clarify the connection of the persistent features with the spin-momentum locking.
Spin magnetic moment of InSb surface polaritons. As shown schematically in fig.2(a), we consider a magnetic field applied alongê x direction (parallel to the surface). The surface plasmon polariton characterized by its conserved in-plane momentum k makes an angle φ withê x (applied field direction) such that φ ∈ [−π, π]. Each such polariton also carries a transverse spin locked to its momentum, depicted in the schematic (spin momentum locking [14]). Figure 2(b) displays the dispersion ω(k ) of polaritons for different angles φ, obtained numerically as described in the methods section. In the absence of magnetic field or with the magnetic field perpendicular to the surface, the surface polaritons are p-polarized and the dispersion is the same as the dispersion for φ = 0 (green curve) for all angles. On the other hand, in presence of magnetic field parallel to the surface, the surface polaritons contain both s-polarized and p-polarized electromagnetic fields, requiring numerical method for calculating the dispersion for arbitrary propagation directions. As shown in the figure 2(b), assuming magnetic field alongê x direction, the polaritons characterized by φ ≥ 0 (with positive spin component along applied magnetic field) are redshifted while those characterized by φ < 0 (with negative spin component along applied magnetic field) are blueshifted. This is further demonstrated in fig. 2(c) where ω(k ) is obtained as a function of angle φ for a fixed k = |k |, for two different values of magnetic field. We numerically find that the energy shift for each k polariton increases linearly with magnetic field in the weak field regime (B 1T) while the dependence is complicated for strong applied fields. All these results strongly indicate that the polaritons have a spin magnetic dipole moment (µ p ) parallel to the transverse spin that interacts with the applied magnetic field. The energy of this Zeeman interaction is described by the Hamiltonian, It follows that the magnetic-field-induced frequency shift (∆ω) for each polariton is: We find that the magnetic moment µ p = |µ p | for each k polariton depends not only on momentum k = |k | but also on the angle φ. This is evident upon closer inspection of fig 2(c) where the energy shift as a function of angle is not sinusoidal but exhibits slight deviation (see maximum and minimum values). In fig. 2(d), the magnetic moment in units of µ B (Bohr magneton) as a function of k is displayed for two different angles φ = −π 4 , 3π 4 showing that µ p (φ) = µ p (φ + π) or µ p (k ) = µ p (−k ). This asymmetry in polaritonic spin magnetic moment and magnetic-field induced energy shift lies at the origin of the asymmetry in the spin angular momentum and heat flux carried by thermally excited ±k polaritons, resulting into PTPS and PPHC parallel to the surface.
We note that the Poynting vector and the spin of polaritonic waves deviate from their usual directions (P k , S ⊥ k ) for large magnetic fields. However, these deviations are small for weak magnetic fields (B 1T) considered above and therefore, this analysis suffices to qualitatively predict the existence of PTPS and PPHC. In the following, we demonstrate them using full fluctuational electrodynamic calculations. PTPS and PPHC in thermal near-field of InSb slab. We compute the spin angular momentum density (Eq.3) and Poynting flux (Eq.2) in thermal near-field of InSb slab in presence of magnetic field of strength 0.5T alongê x direction. Both vacuum and material are assumed to be at thermodynamic equilibrium temperature of T = 300K. Based on the discussion of polaritons in the previous section, we calculate spin-resolved quantities in the sense described schematically in fig.3(a). In particular, the contributions of electromagnetic waves characterized by φ ≥ 0 (k (+) -waves) and those characterized by φ < 0 (k (−) -waves) are calculated separately. Figure 3(b,c,d) demonstrate the frequency spectra of energy density W (ω), spin angular momentum density S(ω) and poynting flux P(ω) at a distance of d = 1µm above the surface of InSb. All these figures depict the separate contributions of k (+) waves (red curves) and k (−) waves (blue curves) along with the sum total (green curves). As evident from fig.3(b), the collective energy density of k + waves is redshifted and that of k − waves is blueshifted similar to polaritons, leading to broadening of total energy density spectrum (green), compared to the spectrum in the absence of magnetic field (black dashed line). The asymmetric overall contributions of k (+) and k (−) waves result into nonzero spin angular momentum density and nonzero Poynting flux at thermal equilibrium i.e. PTPS and PPHC. Note that these persistent quantities contain contributions from not only surface localized polaritons but also other evanescent waves. For instance, another small peak apparent in 3(b) and clearly visible in 3(d) is not related to the polaritons studied in previous section but instead arises from other nonreciprcoal surface waves which make small contribution in comparison to surface plasmon polaritons. Figure 4 describes PTPS spectrum [4(a,b,c)] and also demonstrates that the total frequency-integrated PTPS [4(d)] can compete with the angular momentum density contained in the laser light. For brevity, we focus only on PTPS. As shown in fig 4(a), the electric-type PTPS is evidently much larger than the magnetic-type PTPS. We describe later [ fig 6] that this holds universally for any gyro-electric type nonreciprocal material. Fig 4(b) demonstrates the change in the spectrum as a function of distance from the surface. At each frequency, the sign (direction) of PTPS stays the same while the magnitude decays exponentially as a function of distance from the surface. This also indicates (although not shown separately) that PTPS and PPHC arise from the waves that are evanescent on the vacuum side of the geometry. Fig 4(c) depicts the dependence of the spin angular momentum density on the applied magnetic field. First, the PTPS spectrum broadens as magnetic field is increased from B = 0.5Tê x (light green) to B = 2Tê x (green). Second, while perpendicular magnetic field by itself does not lead to PTPS, it does affect the spectrum observed in presence of parallel magnetic field. This is evident from green and dark green (B = 2T (ê x +ê z )) curves. When magnetic field is at oblique angles to the surface such as this example (dark green), the cyclotron motion in yz-plane due to B x is intercoupled with that in xy-plane due to B z . Due to this intercoupled cyclotron motion of underlying charges, the perpendicular component of magnetic field (which does not lead to PTPS by itself) affects the PTPS spectrum obtained with magnetic field parallel to the surface.
Finally, we plot the total frequency integrated PTPS We note that the total frequency-integrated Poynting fluxP y for this practical example (not shown) is along −ê y direction. The direction of integrated PPHC is related to the underlying cyclotron motion of electrons induced by magnetic field [19]. Inside the bulk of InSb, the cyclotron motions of electrons cancel each other but at the surface, this cancellation is incomplete. The direction of this incomplete cyclotron motion co-incides with the direction of the PPHC.
We further note that the practical example considered here can also lead to unidirectional energy transport because of nonreciprocity [42,46,47].
However, unidirectional transport is not a cause of PTPS and PPHC and this is explained in the following.  fig.3(d) has a predominant contribution along the same direction, both PPHC and PTPS spectra shown in fig. 3(c,d) show that smaller frequencies with bidirectional polaritons also lead to nonzero persistent quantities. While recent works [48] have started to explore the role of nonlocality in context of nonreciprocity and unidirectional transport, we leave nonlocality aside for future work.
Experimental Proposal.
Here we propose an experiment that can provide a direct visual evidence of PPHC using the planar slab considered above. It is well-known that light carrying momentum and spin/orbital angular momentum can exert optical forces and torques on small absorptive particles in its path. Many works have explored this light-matter interaction in optical tweezers [35,49,50] by non-thermal means and in entirely passive systems where forces and torques originate from intrinsic quantum and thermal fluctuations [51][52][53]. We are interested in the latter for non-intrusive (without disturbing thermal equilibrium) detection of the persistent phenomena. We therefore consider a passive system shown schematically in figure 5 where aqueous or fluidic environment covers the gyrotropic InSb thick slab and contains suspended small micrometer size absorptive and nonmagnetic particles. The particles perform Brownian motion about their positions at finite equilibrium temperature. Upon application of magnetic field, the particles experience additional optical force and torque associated with PPHC and PTPS. While the average motional energies of particles remain constant at thermal equilibrium (∼ 1 2 k B T by equipartition law), the additional forces and torques lead to preferential changes in their mean positions and angular orientations which can be detected using a microscope [54,55]. In the following, we estimate these changes with simplifying approximations. A rigorous description of Brownian dynamics [56,57] is beyond the scope of this work and at this point unnecessary since the goal is to merely detect the presence of the persistent thermal photonic phenomena. The particles are absorptive, nonmagnetic (not influenced by presence or absence of magnetic field) and much smaller in size ( µm) and hence dipolar in nature at surface-polariton wavelengths ( 48µm) of InSb. The entire system is at thermal equilibrium room temperature. Magnetic field of 1T is applied alongê x direction resulting into persistent planar heat current P = −P yêy and persistent thermal photon spin S = −S xêx . Our analysis above of planar geometry with infinite transverse extent is valid when both magnetic field and temperature are uniform over an area much larger than polaritonic wavelengths which are of the order of 100µm. Assuming that these conditions are realized over an area of 1cm 2 of InSb slab in an actual experiment which is quite realistic, we extend our fluctuational electrodynamic analysis to calculate the additional forces acting on particles, originating from the local persistent features (not influenced by the edge effects).
The average stochastic optical force on a Brownian particle is written as a sum of the following two terms [58,59]: Here p is the dipole moment of the particle and E is the total electric field at the position of the particle. The first term corresponds to the interaction of fluctuating dipole moment of the particle with thermal fields induced by the particle itself. By calculating induced thermal field using Green's function, it can be shown that this does not lead to any lateral force on the particle because of translational invariance [60]. Note that we are primarily interested in the lateral forces since the perpendicular forces exist in the near-field of all materials [58] and cannot be used to detect PPHC. The second term denotes the interaction of the thermal fields with the induced dipole moment of the particle given by (ω)+ b is the polarizability of the spherical nanoparticle of isotropic permittivity ε = (ω) and radius R, immersed in water of permittivity b ≈ 1.77. This interaction leads to the following simplified expression [35] for the lateral force which occurs only alongê y direction for magnetic field applied alongê x direction: Here, z denotes the distance of the particle from InSb slab surface and c is speed of light. There are many readily available particles such as chalk or milk particles and other commercially available nanoparticles which can be used for the Brownian experiment. For estimation purpose, we consider doped Silicon particles of diameter 1µm, mass density ρ ∼ 2329kg/m 3 and Drude permittivity dispersion ε(ω) = 11.7 − ω 2 p /(ω 2 + iγω) where ω p = 1.3 × 10 14 rad/s and γ = ω p /100. For these parameters, at a distance of z = 2µm from the slab surface, the particles experience linear acceleration of 4µm/s 2 along −ê y i.e. direction of PPHC. While the particles also experience torque due to PTPS, the resulting rotational changes are difficult to observe with the proposed experiment. Nonetheless, the above calculations indicate that there should be a noticeable displacement of Brownian particles along the direction of PPHC which can be viewed using a microscope as depicted in fig. 5. Since the additional lateral movement at thermal equilibrium is not possible in absence of magnetic field or with other homogeneous reciprocal media, its mere presence will be a clear indicator of the persistent thermal photonic phenomena. It can be readily perceived upon seeing through a simple lab microscope as shown in the insets of fig. 5 or by methodically tracking the particle movements.
ix Universal behavior of PTPS and PPHC. We now describe the universal behavior of PTPS and PPHC with generic biansiotropic material types. A bianisotropic medium is often considered in the literature [61,62] to represent a superset of all types of media, more commonly described with following constitutive relations assuming local material response (in the frequency domain): ε, µ are dimensionless permittivity and permeability tensors and ξ, ζ are magneto-electric coupling tensors. Based on the existing literature, we categorize bianisotropic materials into following five well-studied material types: • Isotropic materials: Most naturally existing dielectric or metallic materials with scalar ε, µ and ξ, ζ = 0.
• Uniaxial/biaxial anisotropic materials such as birefringent crystals that have diagonal ε and µ with unequal diagonal entries and ξ, ζ = 0.
Apart from the naturally existing examples given above, there exists a huge range of artificially designed metamaterials with/without bias fields/currents that can effectively provide any combination of these material types [61,62]. It is well-known that such a bianisotropic material is reciprocal [68] when, The material is nonreciprocal if atleast one of these conditions is violated. Table I summarizes PTPS and PPHC for generic bianisotropic material classes with suitable representative examples that describe the presence or absence of PTPS and PPHC. Both isotropic and uniaxial/biaxial anisotropic materials being reciprocal in nature, do not lead to any persistent spin or heat current. The first example considers uniaxial anisotropic material with its anisotropy axis parallel to the surface (breaking the rotational symmetry) and the full calculations confirm the absence of the persistent phenomena. The second and third examples in the table correspond to gyroelectric and gyromagnetic materials having anti-symmetric (nonreciprocal) permittivity ε and permeability µ tensors respectively. For both examples, the gyrotropy axis is assumed to be alongê x which leads to PTPS alongê x and PPHC alongê y direction. It is also found that PTPS of gyroelectric material is mostly electric-type while that of gyromagnetic material is mostly magnetic-type. While the chosen parameters lead to plasmonic enhancement of the persistent phenomena, other parameters (dielectric ε and µ) also show the same (zero or nonzero) features. In the fourth example with gyrotropy axis alongê z (perpendicular to surface), no PPHC is observed. Thê e z components of electric and magnetic type PTPS are nonzero but cancel each other leading to zero total thermal spin. This proves that nonreciprocity is not sufficient to observe PTPS and PPHC.
The fifth and sixth examples consider isotropic, dielectric permittivity and permeability (ε, µ) and diagonal magneto-electric susceptibilities (ξ, ζ). For nonreciprocal susceptibilities (ξ = −ζ T ), it is found that both PPHC and PTPS are zero although electric and magnetic type contributions to PTPS are nonzero. This is qualitatively similar to gyrotropic media with gyrotropy axis perpendicular to the surface. When the off-diagonal components of ξ, ζ are nonzero as considered in the seventh (last) example, PTPS and PPHC parallel to the surface are observed. Interestingly, when y and z-components of E, H fields are coupled, PTPS is alonĝ e y direction while PPHC is alongê x direction. When x and y-components of E, H fields are coupled, PTPS and PPHC are parallel to the surface but they are not necessarily in a specific direction (not tabulated). For a magneto-electric medium, both electric and magnetic contributions to PTPS are comparable to each other as opposed to gyrotropic media where one of them dominates. Figure 6 summarizes the important findings based on this general analysis.
We note that the material parameters ε, µ, ζ, ξ are not entirely arbitrary but follow certain symmetry relations and are also constrained by conditions of causality and passivity [69,70]. The causality constraint leads to Kramer-Kronig relations for frequency dependent parameters and it requires separate examination for different types of materials [71]. The passivity requires x Nonreciprocal materials that lead to PTPS and PPHC that the material matrix, is such that [(−iM )+(−iM ) * T ]/2 is positive definite [69]. We do not discuss the frequency dependence of various parameters here since single frequency calculations are sufficient to describe the nature (existence, directions) of PTPS and PPHC for a given material type. However, we make sure that all the parameters satisfy the constraint of passivity and note that without such constraints, the persistent phenomena can be incorrectly deduced for reciprocal systems that are non-passive (nonequilibrium).
The seven examples above and the analysis presented here are sufficient to predict the presence or absence of the PTPS and PPHC and their nature (directions, electric/magnetic type PTPS) for any practical example of a bianisotropic material.

III. CONCLUSION
Modern thermal photonics utilizes fluctuational electrodynamic paradigm to explore new phenomena (near-field radiative heat transfer [40]) and new regimes (nonreciprocity [20], nonlinearities [72] and nonequilibrium [73]) which are inaccessible to older paradigm of radiometry and Kirchhoff's laws. And yet, thermal spin photonics is so far limited to inquiries xi based on Kirchhoff's laws [7][8][9]. This work demonstrates intriguing spin angular momentum related thermal radiation phenomena in the near-field of nonreciprocal materials analyzed within fluctuational electrodynamic paradigm.
It paves the way for new fundamental and technological avenues in thermal spin photonics. In particular, it will be useful in the near future for shaping spin-angular-momentum related radiative heat transport phenomena such as our recent work on circularly polarized thermal light sources [11].
Our work revealed that the spin-momentum locking of thermally excited evanescent waves plays a fundamental role in facilitating the surprising thermal equilibrium features of PTPS and PPHC. The connection between the spin-momentum locking and the radiative heat transfer is important for exploring new ways of achieving directional heat transport at the nanoscale. We found that the surface polaritons of gyrotropic materials can carry spin magnetic moment which invites separate related studies of spin-dependent quantum plasmonics and spin quantization. We proposed an experiment based on Brownian motion that can provide a visual confirmation of the persistent phenomena. Currently, there are no experiments probing such intriguing nonreciprocal thermal fluctuations and heat transport effects. Also, the experimental detection of the predicted surprising effects is important from the perspective of thermodynamic revalidation of fundamental understanding of nonreciprocal systems.
We described the universal behavior of the thermal spin photonic phenomena with a comprehensive analysis of key classes of nonreciprocal materials namely, gyroelectric, gyromagnetic and magneto-electric media. This general analysis motivates similar studies of thermal radiation [74], radiative heat transfer [40], Casimir forces/torques [39] from generic bianisotropic materials including largely unexplored material types in this context such as topological insulators, multiferroic and magentoelectric materials.
The theoretical framework and the Green's function produced here can also be used to study environment-assisted quantum nanophotonic phenomena such as Forster resonance energy transfer [26], atomic transition shifts [27] with general, bianisotropic materials. We leave all these promising directions of research aside for future work.
Acknowledgments This work was supported by the U.S. Department of Energy, Office of Basic Energy Science under award number DE-SC0017717, DARPA Nascent Light-Matter Interaction program and the Lillian Gilbreth Postdoctoral Fellowship program at Purdue University (C.K.).

IV. METHODS
Derivation of Green's function. The Green's function relating vector potential at r 1 = (R 1 , z 1 ) to source current at r 2 = (R 2 , z 2 ) in vacuum where R = (x, y) denotes planar co-ordinates, is: k = (k , ±k z ) is the total wavevector consisting of conserved parallel component k and perpendicular z-component in vacuum ±k z . The (+) and (−) signs denote waves going away from and towards the interface respectively. It follows from Maxwell's equations that they satisfy the dispersion relation k 2 +k 2 z = k 2 0 = (ω/c) 2 where k = |k | is real and k z can be real (k < k 0 ) or complex valued (k > k 0 ). For simplicity, we write k = (k cos φ, k sin φ) where φ is the angle subtended by k with x-axis. Assuming z 1 > z 2 , the integrand g(k , z 1 , z 2 ) tensor is written using the s, p-polarization vectors (ê s ,ê p ): The polarization vectorsê j± for j = s, p with ± denoting waves going along ±ê z directions are: The Fresnel reflection coefficient r jk for j, k = [s, p] describes the amplitude ofê j -polarized reflected light due to unit amplitudeê k -polarized incident light. The Green's function above consists of two parts corresponding to the trajectories of electromagnetic waves generated at the source position r 2 and arriving at r 1 either directly (g 0 ) or upon reflection from the interface (g ref ). The Green's function in Eq. 15 xii is derived for z 1 ≥ z 2 . For z 1 < z 2 , only the vacuum part is modified to g 0 = e −ikz(z1−z2) [ê s−ê T s− + e p−ê T p− ]. The Fresnel reflection coefficients can be obtained experimentally or theoretically.
Fresnel reflection coefficients and polaritonic dispersion. We develop a tool to compute Fresnel reflection coefficients for a generic, homogeneous medium that can be described using the following constitutive relations assuming local response (in the frequency domain): ε, µ are dimensionless permittivity and permeability tensors and ξ, ζ are magneto-electric coupling tensors. For isotropic materials, ξ, ζ = 0 and ε, µ are scalars. For gyro-electric (magneto-optic) and gyro-magnetic media, the tensors ε and µ respectively have off-diagonal components and ξ, ζ = 0. The tensors ξ, ζ are nonzero for magneto-electric media. By writing electromagnetic fields inside the material as [E, µ0 ε0 H] T e i(k ·R+kzz−iωt) and using above constitutive relations in Maxwell's equations, we obtain the following dimensionless dispersion equation for waves inside the material: Here, 6 × 6 material tensor M describes the constitutive relations and M k corresponds to the curl operator acting on plane waves. Because of the generality of this problem, we obtain k z numerically by solving det(M + M k (k z )) = 0 for given (k , φ). Depending on the nature of the material, there can be two (for isotropic media) or four (for anisotropic media) solutions of k z corresponding toê z -propagation of electromagnetic waves. Overall, there are four eigensolutions spanning the null-space of the matrix M + M k (k z ), out of which two solutions correspond to waves propagating in −ê z direction (transmitted waves in our geometry). The four Fresnel reflection coefficients are then obtained by matching the tangential components at the interface (E x , E y , H x , H y ) of incident and reflected fields with the transmitted fields. Here, the transmitted fields are written in the basis of former two null-space solutions while the incident and reflected fields are written in the basis ofê s ,ê p -polarizations (Eq. 16). This procedure is also extended in this work to compute the polaritonic dispersion (ω(k )) of surface polaritons that decay on both sides of the interface. While that calculation does not involve Fresnel coefficients, the boundary conditions again lead to a homogeneous, linear problem of the form M p (ω, k )X = 0 where X contains the coefficients describing the decomposition of polaritonic fields into four eigenstates (s, p-polarizations in vacuum and two −ê z -propagating solutions inside the medium) at the interface. By numerically solving det(M p (ω, k )) = 0, the polaritonic dispersion ω(k ) is obtained. The associated null-space describes the polaritonic fields. Note that since k is assumed to be real-valued and non-decaying, ω is complex-valued with the imaginary part describing the finite lifetime (quality factor) of the polaritons. We consider a semi-infinite half-space of a generic bianisotropic medium at thermal equilibrium with vacuum. To analyze the thermal radiation on the vacuum side of the geometry, we derive the Green's function, equilibrium correlations of vector potential (fluctuation dissipation relation) and equilibrium correlations of electromagentic fields. Finally, we provide semi-analytic expressions for spin angular momentum density and Poynting flux perpendicular and parallel to the surface.

I. DERIVATION OF GREEN'S FUNCTION
The vector potential A(r 1 ) at r 1 produced by source current density J(r 2 ) located at r 2 can be calculated using Green's function G(r 1 , r 2 ) using the relation: Since we use Landau gauge E = iωA, this is same as the electric-type Green's function that is commonly employed in the literature in the form of following equation: where p is polarization (dipole moment) density. Since the derivation of this Green's function is quite well-known for isotropic media, we do not reproduce that derivation here but focus mainly on its extension to the case of bianisotropic half-space considered in the manuscript. We use Weyl's angular spectrum representation of Green's function. We write the position vectors as r j = (R j , z j ) with transverse co-ordinates R j = (x j , y j ) for j = [1, 2] and wave-vectors k = (k , k z ) with transverse wavevector k = k (cos φê x + sin φê y ) where we introduce φ for simplicity of expressions below. In vacuum, the dispersion relation k 2 + k 2 z = k 2 0 = (ω/c) 2 follows from Maxwell's equation, where k is real-valued and k z can be real (for k ≤ k 0 ) or complex-valued (for k > k 0 ). In vacuum, the electric field at r 1 produced by the dipole moment p = dδ(r − r 2 ) located at r 2 is written in the angular spectrum representation for z 1 ≥ z 2 as: The polarization vectorsê j± for j = s, p with ± denoting waves going along ±ê z directions are: For z 1 < z 2 , the integrand is modified and contains the term e −ikz(z1−z2) [ê s− (ê s− · d) +ê p− (ê p− · d)]. For consistency, we will stick with z 1 ≥ z 2 in the following discussion. The scattered/reflected field is calculated by considering the reflection of the incident field at the interface (z = 0). Since the waves propagate in −ê z direction to reach the interface, the incident field will be of the form e ikzz2 [ê s− (ê s− · d) +ê p− (ê p− · d)]. It undergoes reflection at the interface where polarization vectors change toê s− → r ssês+ + r psêp+ andê p− → r spês+ + r ppêp+ . The Fresnel reflection coefficient r jk for j, k = [s, p] describes the amplitude ofê j -polarized reflected light due to unit amplitudê e k -polarized incident light. For isotropic media, cross-polarization Fresnel coefficients r sp , r ps are zero which simplifies the calculation. But they are not necessarily zero for general bianisotropic media. The reflected field then acquires ii an additional phase of e ikzz1 upon reaching the position r 1 along with the overall transverse phase accrual same as e ik ·(R1−R2) . This results in the scattered/reflected field at r 1 given below: E ref (R 1 , z 1 ) = ω 2 µ 0 d 2 k (2π) 2 e ik ·(R1−R2) i 2k z e ikz(z1+z2) [(r ssês+ + r psêp+ )(ê s− · d) + (r spês+ + r ppêp+ )(ê p− · d)] By writing the total field E(r 1 ) = E 0 (r 1 ) + E ref (r 2 ) and using Eq 2, the Green's function G(r 1 , r 2 ) is derived for the geometry considered in the manuscript.

II. DERIVATION OF FLUCTUATION DISSIPATION THEOREM (FDT)
We follow Landau's discussion in Ref. [1] (Statistical Physics Part 2, Chapter 8) to obtain the vector potential correlations in the vacuum half-space when both the vacuum and the material half-spaces are at the same thermodynamic temperature T (FDT of first kind). Let's consider the linear response theory developed by Kubo. In this theory, we consider a discrete set of quantities denoted by x a for (a = 1, 2, ...) which describe the behavior of the system under certain external interactions. These interactions are described by external forces f a such that interaction energy has the form: The quantities x a are further related to to the forces f a through linear generalized susceptibilities α ab (r, r ) (linear response). In the Fourier domain, they can be written as: x a (ω, r) = b α ab (ω; r, r )f b (ω, r )d 3 r The spectral distribution of the fluctuating quantities x a (ω, r) is related to the generalized susceptibilities by Kubo's fluctuation dissipation relation given by: x a (r, ω)x * b (r , ω ) = α ab (ω; r, r ) − α * ba (ω; r , r) 2i For electromagnetic fields, x a (ω, r) → A j (ω, r) (j = x, y, z component of vector potential). The interaction with the externally induced current is given by V = −j · A where j(ω, r) is the generalized force. Since vector potential and current density are related by the Green's function: