Angular momentum radiation from current-carrying molecular junctions
AAngular momentum radiation from current-carrying molecular junctions
Zu-Quan Zhang, ∗ Jing-Tao L¨u, † and Jian-Sheng Wang ‡ Department of Physics, National University of Singapore, Singapore 117551, Republic of Singapore School of Physics and Wuhan National High Magnetic Field Center,Huazhong University of Science and Technology, 430074 Wuhan, P. R. China (Dated: 15 November 2019)We consider the radiation of angular momentum (AM) from current-carrying molecular junctions.Using the nonequilibrium Green’s function method, we derive a convenient formula for the AMradiation and apply it to a prototypical benzene molecule junction. We discuss the selection rules forinelastic transitions between the molecular angular momentum eigenstates due to a 6-fold rotationalsymmetry. Our study provides important insights into the generation of light with AM from DC-biased molecular junctions.
Introduction.–
The atomic scale interaction of nonequi-librium electrons with light is the key to develop elec-trically driven single molecular light sources for sensing,spectroscopy and chemical reactions [1–3]. The high spa-tial resolution and local field enhancement offered bythe tip of a scanning tunneling microscope provide anideal platform to investigate this problem [4–9]. Recentyears have witnessed tremendous progress in this direc-tion. By analyzing the light emission spectra, a varietyof physical and chemical information can be deduced,including vibrational coupling [7, 10], coherent inter-molecular dipole interaction [8, 11, 12], plasmon-excitoncoupling [13, 14], anti- and super-bunching photon statis-tics [15, 16], charge and spin state emission [17, 18]. The-oretical understanding of these effects relies on methodsdeveloped in quantum optics and quantum transport [18–27].The coupling of electron orbital motion with its spinleads to spin-orbit interaction, which is of vital impor-tance in spintronics [28], topological physics [29, 30] andso on. Spin-orbit coupling is also responsible for the chi-ral induced spin selectivity in electron transport throughmolecules [31, 32]. However, the effect of electron or-bital motion on single molecule electro-luminescence is,to a large extent, unexplored. In this work, based onnonequilibrium Green’s function (NEGF) method [33–36], we develop a microscopic theory to study electri-cally driven angular momentum emission from a singlemolecule. We consider a prototypical benzene molecule,which has well-defined orbital angular momentum statesin isolated situation. We illustrate the underlying mecha-nism as inelastic electronic transition between states withdifferent orbital AM. This is in contrast to the opticalapproach by passing normal light through constructedoptical structures [37, 38].
Theory.–
To consider an open system for light emis-sion, we decompose it into four parts: a molecular systemas a quantum emitter, the coupling of the quantum emit-ter with the radiation field, the radiation field itself, theleads and their couplings with the quantum emitter forpumping energy and electrons into the quantum emitter.
FIG. 1. Light emission from a benzene molecule junction.The metal leads L and R are connected to two carbon atomsin the ortho position. Meta and para positions correspond tolead R connecting to 3 and 4, respectively. We use a tight-binding (TB) model combined with thePeierls substitution [39] to describe the central moleculeand its coupling with the radiation field, written as H t = (cid:88) (cid:104) ij (cid:105) c † i t ij c j e iθ ij , (1)where (cid:104) ij (cid:105) denotes the nearest-neighbor (NN) sites i and j , t ij is the NN hopping parameter, c † i ( c i ) is the electroncreation (annihilation) operator at site i . The phase θ ij = e (cid:126) (cid:82) r i r j A · d l represents the coupling to the radiation field.Here, e ≈ − . × − C is the electron’s charge, A isthe vector potential, r i and r j are the positions of sites i and j respectively. Expanding e iθ ij in terms of A upto first order, we can write Eq. (1) into two terms H t ≈ H t + H int , with H t for the noninteracting electrons and H int for the coupling of the electrons with the radiationfield. They are given by H t = (cid:88) (cid:104) ij (cid:105) c † i t ij c j , (2) H int ≈ (cid:88) (cid:104) ij (cid:105) (cid:88) k (cid:88) µ = x,y,z M kµij c † i c j A µ ( r k ) , (3)where M kµij = i e (cid:126) t ij ( r i − r j ) µ ( δ ki + δ kj ) is the electron-photon coupling matrix element. We use Greek letters to a r X i v : . [ c ond - m a t . m e s - h a ll ] N ov represent components of the Cartesian coordinates, i.e., µ = x, y, z .Hamiltonian of the radiation field is H rad = 12 (cid:90) d r (cid:0) ε E ⊥ + 1 µ B (cid:1) , (4)where ε and µ are the vacuum permittivity and perme-ability, respectively. Here, we adopt the Coulomb gaugewith ∇ · A = 0, thus the transverse electric field is givenby E ⊥ = − ∂ t A , and the magnetic field is B = ∇ × A .We restrict our discussion to the far-field radiation hereand neglect the longitudinal electric field, which is im-portant only in the near-field heat transfer [40–42]. Theeffect of the leads and their couplings with the moleculeare included by the self-energies, as shown below.The energy and AM flux of electromagnetic field canbe obtained from [43, 44] S = 1 µ (cid:10) : E ⊥ × B : (cid:11) , (5) ←→ M = (cid:104) : r × ←→ T : (cid:105) , (6)where (cid:104) : AB : (cid:105) denotes normal order of operators AB when taking the ensemble average, which removes thezero-point motion contribution, and ←→ T is the Maxwellstress tensor with T µν = δ µν ( ε E + µ − B ) − ε E µ E ν − µ − B µ B ν . Equations (5) and (6) can be expressed in terms of photon Green’s function (GF) [45] S µ ( r ) = (cid:15) µνγ (cid:15) γδξ µ (cid:90) + ∞ dω π (cid:126) ω × Re (cid:20) − ∂∂x (cid:48) δ D <νξ ( r , r (cid:48) , ω ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) r (cid:48) → r . (7)Einstein summation rule is used here, and (cid:15) µνγ is theLevi-Civita symbol, D <νξ ( r , r (cid:48) , ω ) is the photon’s lesserGF in the frequency domain. Relevant quantities inEq. (6) are written as (cid:104) : E µ E ν : (cid:105) = 2 i (cid:126) (cid:90) ∞ dω π ( (cid:126) ω ) Re (cid:2) D <µν ( r , r , ω ) (cid:3) , (8a) (cid:104) : B µ B ν : (cid:105) = i (cid:126) (cid:90) ∞ dω π (cid:15) µγξ (cid:15) νγ (cid:48) ξ (cid:48) × Re (cid:20) ∂∂x γ ∂∂x (cid:48) γ (cid:48) D <ξξ (cid:48) ( r , r (cid:48) , ω ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) r (cid:48) → r . (8b)The GFs are obtained following the standard NEGFformalism. The retarded GF is solved by the Dyson equa-tion D r = d r + d r Π r D r , and d r is the free space photonGF [46]. The lesser GFs are obtained by the Keldyshequation D < = D r Π < D a , with D a = ( D r ) † . We use therandom phase approximation to calculate the interactingself-energyΠ <µν ( r i , r j , ω ) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) Tr (cid:104) M iµ g < ( E ) M jν g > ( E − (cid:126) ω ) (cid:105) , (9)where Tr[ · · · ] means trace over the electron degrees offreedom, g < ( > ) = g r Σ < ( > )leads g a is the lesser (greater) GFsfor the noninteracting electrons. Here, Σ < ( > )leads is the lesser(greater) self-energy due to electron’s coupling with theleads. The leads are in their respective equilibrium statesand the self-energies follow the fluctuation-dissipationtheorem, for lead β , Σ <β = if β Γ β and Σ >β = i ( − f β )Γ β .Here, f β ( E, µ β ) = 1 (cid:46)(cid:104) exp (cid:0) E − µ β k B T β (cid:1) + 1 (cid:105) is the Fermi dis-tribution function, µ β is the chemical potential, T β isthe temperature, k B is the Boltzmann constant, andΓ β = − (cid:2) Σ rβ (cid:3) is the coupling-weighted spectrum ofthe lead.To calculate the total energy and AM radiation, wechoose a large spherical surface enclosing the molecule,and perform the surface integral P = (cid:73) S · d A , (10) d L dt = (cid:73) ←→ M · d A , (11) where d A = ˆ R dA , with ˆ R = R /R denoting the unitnormal vector of the spherical surface dA with radius R . In the far-field region (with R much larger than thephoton wavelength λ and the central molecule’s size ∼ a ),we get simplified expressions for Eq. (10) and Eq. (11) as P = − (cid:90) ∞ dω π (cid:126) ω πε c Im (cid:2) Π tot ,<µµ ( ω ) (cid:3) , (12) dL γ dt = (cid:90) ∞ dω π (cid:126) ω πε c (cid:15) γµν Re (cid:2) Π tot ,<µν ( ω ) (cid:3) , (13)with Π tot ,<µν ( ω ) = (cid:80) ij Π <µν ( r i , r j , ω ). Einstein summa-tion rule is used here. We have performed the solid angleintegration within the monopole approximation, i.e., ne-glecting the molecular size considering ( a/R ) (cid:28)
1. Weobserve that the emission power is related to the trace,while AM emission rate to the antisymmetric part of thetensor, Π tot ,< . Eq. (13) is the main result of this pa-per. It serves as a simple formula for calculating far-fieldradiation of optical AM from biased molecular systems. Application.–
The theory we develop is quite general.We now apply it to a prototypical benzene molecule junc-
FIG. 2. (a) Selection rules for the emission of optical an-gular momentum between energy levels of the TB benzenemolecule. (b) Light emission with leads coupled directly totwo eigenmodes of the benzene molecule in high bias regime.¯Γ is the weak lead coupling to the modes. (c) Orbital resolvedelectron density of states (DOS) of the benzene molecule whencoupled to the leads in ortho position shown in Fig. 1, withΓ = 0 . − Im (cid:2)(cid:101) g rll ( E ) (cid:3) /π for energylevel l . tion shown in Fig. 1. We take the NN hopping parameteras − t ij = t = 2 . a = 1 . L = Γ R = Γ.The isolated benzene molecule has a C rotationalsymmetry. We put the molecule in the x - y plane withsite positions x j = a cos(2 πj/ y j = a sin(2 πj/ j = 1 , , . . . ,
6. The orbital eigen energies are E l = − t cos(2 πl/ l = 0 , ± , ± ,
3. We have neglectedthe spin degeneracy here. When coupled to the twometal leads, the energy levels are broadened, but the de-generacy is not lifted. The orbital resolved density ofstates (DOS) is shown in Fig. 2(c). Here the mode spaceGF (cid:101) g r is related to the real space GF via the unitarytransformation (cid:101) g r = U † g r U , with U jm = e i πjm/ / √ j, m = 1 , , . . . ,
6. Note that the mode index m is uniqueonly modulo 6, thus 6 is the same as 0, and 5 is the sameas −
1. Also, we use the notation that an operator de-noted as O in real space is written as (cid:101) O in mode space,with (cid:101) O = U † OU .The mechanism of AM emission can be understood byconsidering the selection rules in mode space. Definingthe electron velocity matrix v µij = t ij ( r i − r j ) µ / (cid:126) , wehave v µ = ie (cid:80) k M kµ , which has the C rotational sym-metry. For the coordinate system we choose here, the C symmetry leads to the relations (cid:101) v xnm (cid:101) v xmn = (cid:101) v ynm (cid:101) v ymn and (cid:101) v xnm (cid:101) v ymn = i ∆ mn (cid:101) v xnm (cid:101) v xmn . Here, ∆ mn = sgn( m − n ) if | m − n | = 1, and ∆ = − ∆ = 1, otherwise ∆ mn = 0.We consider the simple case where the leads couple re-spectively to only two eigenmodes of the molecule in thehigh bias regime µ L (cid:28) E n < E m (cid:28) µ R [Fig. 2(b)].Due to the high bias setting, we have (cid:101) g 0. In the limit ¯Γ → 0, we get fromEq. (12) P = − ω mn e ( (cid:101) v xnm (cid:101) v xmn + (cid:101) v ynm (cid:101) v ymn ) / (3 πε c ),and dL z /dt = iω mn e ( (cid:101) v xnm (cid:101) v ymn − (cid:101) v ynm (cid:101) v xmn ) / (3 πε c ) fromEq. (13), with (cid:126) ω mn = E m − E n [45]. Using the relationsof the velocity matrix due to the C symmetry, we get dL z /dtP = ∆ mn ω mn . This result is reminiscent of the Eq. (19)of a recent work for the classical case of AM radiationfrom a single electron performing circular motion witha constant frequency [47]. Since every emitted photoncarries an energy (cid:126) ω mn , the number of photons emittedper unit time is dN/dt = P/ ( (cid:126) ω mn ). Thus, the AM peremitted photon is ∆ L = dL z /dtdN/dt = ∆ mn (cid:126) , which is theselection rules shown in Fig. 2(a).The light emission for real space coupling in Fig. 1is a combination of all the possible processes shown inFig. 2(a). This is tuned by the applied bias. Significantlight emission between two energy levels of the moleculeis possible when the energy levels enter into the bias win-dow, restricted by the selection rules. Since the degen-erate energy levels ( l = ± , ± 2) are broadened but notsplit when the molecule couples to the leads, they will en-ter into or out of the bias window simultaneously. Lightemission from inelastic transition l = 2 → l = 1 is ac-companied by emission from transition l = − → l = − µ L and µ R in the ortho position. We observe a four-line-segmentfeature (FLSF) where the AM radiation is large when onechemical potential is in resonance with the eigenstates at ± t and the bias window covers the energy range [ − t, t ].The AM radiation is quite small in other regions. Toanalyze this resonant effect, we show in Fig. 3(b) thefrequency/energy resolved spectrum of the AM radiation(Eq. (13) without frequency integration) for the resonantcase. There is a large peak at (cid:126) ω = 2 t , which impliesthat the FLSF in Fig. 3(a) is contributed mainly fromthe radiative transition processes l = 2 → l = 1 and l = − → l = − E = ± t . We find thatthe cross correlations between the degenerate states, suchas (cid:101) g r ( E ) and (cid:101) g r ( E ), are important in resulting the netAM radiation than the correlations between nondegener-ate states, such as (cid:101) g r ( E ) and (cid:101) g r ( E ). We set the latterto 0 for simplicity. With these simplifications, for theresonant peak at µ L = − t , we get from Eq. (13) in the . . . . . . . . . . − . − . − . − . − . − . − . . . . − . − . − . − . − . − − µ L /t − − µ R / t − . − . . . . (a) A M r a d i a t i o n s p ec tr u m ¯ hω/t (b) ( d L z / d t ) / J µ L /t (i)(ii)(iii)(iv) (c) ( d L z / d t ) / J µ L /t orthometapara (d) FIG. 3. (a) Intensity plot of angular momentum radiationrate normalized by ( dL z /dt ) /J as a function of chemical po-tentials µ L /t and µ R /t . (b) Frequency resolved spectrum ofthe angular momentum radiation, normalized relative to themaximum value at (cid:126) ω = 2 t , for µ L = − t and µ R = 6 eV. (c)Line cut of the plot in (a) at µ R = 4 eV. (i)(ii) are results cal-culated from Eq. (13) at zero temperature, while (iii)(iv) areresults from Eq. (14). Γ = 0 . . µ R = 4 eV. For (a)(b)(d),Γ = 0 . T = 300 K. zero temperature limit [45] dL z dt ≈ J (Γ / ( µ L + t ) + (Γ / , (14)with J = √ π tα ( v /c ) , v = at/ (cid:126) , the fine structureconstant α = e / (4 πε (cid:126) c ), and the light speed c . Theapproximate expression of Eq. (14) agrees well with nu-merical results from Eq. (13) [Fig. 3(c)]. Eq. (14) impliesthat the height of the resonance peak is a constant atzero temperature, and its width is characterized by Γ / Conclusion.– In summary, using the NEGF method,we have developed a theoretical framework to study AMradiation from current-carrying molecular junctions. Asan application, the theory identifies from a quantumviewpoint that electrons of a ring-like benzene moleculeemit light with AM due to radiative transitions betweendifferent angular momentum states. Due to asymmet-ric couplings to the leads and electron tunneling betweendegenerate energy states with opposite angular momen-tum, large resonant effect with the bias potential wasdiscovered. Our theory can be straightforwardly appliedto more realistic chiral molecules.Z.Q.Z. and J.S.W. acknowledge the support by MOEtier 2 Grant No. R-144-000-411-112 and FRC GrantNo. R-144-000-402-114. J.T.L. is supported by the Na-tional Natural Science Foundation of China (Grant No.21873033). ∗ [email protected] † [email protected] ‡ [email protected][1] M. Galperin and A. Nitzan, Phys. Chem. Chem. Phys. , 9421 (2012).[2] S. V. Aradhya and L. Venkataraman, Nat. Nanotechnol. , 399 (2013).[3] K. Kuhnke, C. Große, P. Merino, and K. Kern, Chem.Rev. , 5174 (2017).[4] R. Berndt, J. K. Gimzewski, and P. Johansson, Phys.Rev. Lett. , 3796 (1991).[5] R. Berndt, R. Gaisch, J. K. Gimzewski, B. Reihl, R. R.Schlittler, W. D. Schneider, and M. Tschudy, Science , 1425 (1993).[6] J. Aizpurua, G. Hoffmann, S. P. Apell, and R. Berndt,Phys. Rev. Lett. , 156803 (2002).[7] X. H. Qiu, G. V. Nazin, and W. Ho, Science , 542(2003).[8] Z.-C. Dong, X.-L. Guo, A. S. Trifonov, P. S. Dorozhkin,K. Miki, K. Kimura, S. Yokoyama, and S. Mashiko,Phys. Rev. Lett. , 086801 (2004).[9] N. L. Schneider, G. Schull, and R. Berndt, Phys. Rev.Lett. , 026601 (2010).[10] C. Chen, P. Chu, C. A. Bobisch, D. L. Mills, and W. Ho,Phys. Rev. Lett. , 217402 (2010).[11] Y. Zhang, Y. Luo, Y. Zhang, Y.-J. Yu, Y.-M. Kuang,L. Zhang, Q.-S. Meng, Y. Luo, J.-L. Yang, Z.-C. Dong,and J. G. Hou, Nature , 623 (2016).[12] H. Imada, K. Miwa, M. Imai-Imada, S. Kawahara,K. Kimura, and Y. Kim, Nature , 364 (2016).[13] Y. Zhang, Q.-S. Meng, L. Zhang, Y. Luo, Y.-J. Yu,B. Yang, Y. Zhang, R. Esteban, J. Aizpurua, Y. Luo,J.-L. Yang, Z.-C. Dong, and J. G. Hou, Nat. Commun. , 15225 (2017).[14] H. Imada, K. Miwa, M. Imai-Imada, S. Kawahara,K. Kimura, and Y. Kim, Phys. Rev. Lett. , 013901(2017).[15] L. Zhang, Y.-J. Yu, L.-G. Chen, Y. Luo, B. Yang, F.- F. Kong, G. Chen, Y. Zhang, Q. Zhang, Y. Luo, J.-L.Yang, Z.-C. Dong, and J. G. Hou, Nat. Commun. , 580(2017).[16] C. C. Leon, A. Rosawska, A. Grewal, O. Gunnarsson,K. Kuhnke, and K. Kern, Sci. Adv. , eaav4986 (2019).[17] B. Doppagne, M. C. Chong, H. Bulou, A. Boeglin,F. Scheurer, and G. Schull, Science , 251 (2018).[18] K. Miwa, H. Imada, M. Imai-Imada, K. Kimura,M. Galperin, and Y. Kim, Nano Lett. , 2803 (2019).[19] M. Galperin and A. Nitzan, Phys. Rev. Lett. , 206802(2005).[20] N. L. Schneider, J. T. L¨u, M. Brandbyge, and R. Berndt,Phys. Rev. Lett. , 186601 (2012).[21] J.-T. L¨u, R. B. Christensen, and M. Brandbyge, Phys.Rev. B , 045413 (2013).[22] F. Xu, C. Holmqvist, and W. Belzig, Phys. Rev. Lett. , 066801 (2014).[23] K. Kaasbjerg and A. Nitzan, Phys. Rev. Lett. ,126803 (2015).[24] L.-L. Nian, Y. Wang, and J.-T. L¨u, Nano Lett. , 6826(2018).[25] J. Kr¨oger, B. Doppagne, F. Scheurer, and G. Schull,Nano Lett. , 3407 (2018).[26] M. Parzefall and L. Novotny, Rep. Prog. Phys. , 112401(2019).[27] S. Mukamel and M. Galperin, arXiv preprintarXiv:1909.04829 (2019).[28] D. Awschalom and N. Samarth, Physics , 50 (2009).[29] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010).[30] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011).[31] K. Ray, S. P. Ananthavel, D. H. Waldeck, and R. Naa-man, Science , 814 (1999).[32] S. Dalum and P. Hedeg˚ard, Nano Lett. , 5253 (2019).[33] H. Haug and A.-P. Jauho, Quantum Kinetics in Trans-port and Optics of Semiconductors (Springer, Berlin,1996). [34] H. Bruus and K. Flensberg, Many-body quantum theoryin condensed matter physics: an introduction (OxfordUniversity Press, 2004).[35] J.-S. Wang, J. Wang, and J. T. L¨u, Eur. Phys. J. B ,381 (2008).[36] J.-S. Wang, B. K. Agarwalla, H. Li, and J. Thingna,Front. Phys. , 673 (2014).[37] S. Zhang, H. Wei, K. Bao, U. H˚akanson, N. J. Halas,P. Nordlander, and H. Xu, Phys. Rev. Lett. , 096801(2011).[38] Y. Gorodetski, A. Drezet, C. Genet, and T. W. Ebbesen,Phys. Rev. Lett. , 203906 (2013).[39] M. Graf and P. Vogl, Phys. Rev. B , 4940 (1995).[40] A. I. Volokitin and B. N. J. Persson, Rev. Mod. Phys. , 1291 (2007).[41] Z.-Q. Zhang, J.-T. L¨u, and J.-S. Wang, Phys. Rev. B , 195450 (2018).[42] J.-S. Wang, Z.-Q. Zhang, and J.-T. L¨u, Phys. Rev. E , 012118 (2018).[43] S. M. Barnett, J. Opt. B: Quantum Semiclass. Opt. ,S7 (2002).[44] M. Janowicz, D. Reddig, and M. Holthaus, Phys. Rev.A , 043823 (2003).[45] See Supplemental Material for details.[46] O. Keller, Quantum theory of near-field electrodynamics (Springer, Berlin, Germany, 2012).[47] M. Katoh, M. Fujimoto, H. Kawaguchi, K. Tsuchiya,K. Ohmi, T. Kaneyasu, Y. Taira, M. Hosaka, A. Mochi-hashi, and Y. Takashima, Phys. Rev. Lett. , 094801(2017).[48] M. F. Maghrebi, A. V. Gorshkov, and J. D. Sau, Phys.Rev. Lett. , 055901 (2019).[49] C. Khandekar and Z. Jacob, Phys. Rev. Appl. , 014053(2019).[50] A. Ott, P. Ben-Abdallah, and S.-A. Biehs, Phys. Rev. B , 205414 (2018). Supplemental Material: Angular momentum radiation from current-carryingmolecular junctions Zu-Quan Zhang, Jing-Tao L¨u and Jian-Sheng Wang DERIVATION OF THE RADIATION FORMULAS In this section, we give some details of the derivations of Eq. (12) and Eq. (13) for energy radiation and AMradiation respectively in the main text. The photon GF is defined as D µν ( r , r (cid:48) ; τ, τ (cid:48) ) = − i (cid:126) (cid:104) T τ A µ ( r , τ ) A ν ( r (cid:48) , τ (cid:48) ) (cid:105) ,where the operators are in the Heisenberg representation, and T τ is the time-order operator on the Keldysh contour.The required lesser GF in Eqs. (7)-(8b) in the main text is obtained by the Keldysh equation D < = D r Π < D a , with D a = ( D r ) † . Specifically, it is given by D < ( r , r (cid:48) , ω ) = d r ( r , r i , ω ) χ < ( r i , r j , ω ) d a ( r j , r (cid:48) , ω ) , (S1)with χ < = ε − Π < ( ε † ) − . Here, the matrix multiplication in coordinate component subscripts ( µ, ν = x, y, z ) andsummation over sites i, j are implied. We have defined the matrix function ε = − Π r d r and ε † = − d a Π a , with for the identity matrix. Excitations of the photons and screening effect are included in χ < ( r i , r j , ω ). The retardedphoton GF in free space is d rµν ( r, ω ) = − e qr πε rc (cid:0) δ µν − x µ x ν r (cid:1) − πε rc (cid:16) e qr − r q − e qr rq (cid:17)(cid:0) δ µν − x µ x ν r (cid:1) , (S2)where c is the speed of light in free space, and q = i ωc .In calculating the energy radiation and AM radiation, the normal order is introduced to remove the contributionof zero-point motion. The vector field operator is split into a positive frequency part and a negative frequency part, A µ = A (+) µ + A ( − ) µ . The rule of normal order is to put operators with positive frequency part on the right side ofnegative frequency part, i.e., (cid:10) : A (+) µ A ( − ) µ : (cid:11) = (cid:10) A ( − ) µ A (+) µ (cid:11) . We can write − i (cid:126) (cid:104) : A µ ( r t ) A ν ( r (cid:48) t (cid:48) ) : (cid:105) = D >µν ( r , r (cid:48) ; t, t (cid:48) ) + D ( N ) νµ ( r (cid:48) , r ; t (cid:48) , t ) − D ( A ) µν ( r , r (cid:48) ; t, t (cid:48) ) . (S3)Here the normal order and anti-normal order correlation functions are written as D ( N ) µν ( r , r (cid:48) ; t, t (cid:48) ) = − i (cid:126) (cid:10) A ( − ) µ ( r t ) A ν ( r (cid:48) t (cid:48) ) (+) (cid:11) , D ( A ) µν ( r , r (cid:48) ; t, t (cid:48) ) = − i (cid:126) (cid:10) A (+) µ ( r t ) A ν ( r (cid:48) t (cid:48) ) ( − ) (cid:11) . (S4)Applying E ⊥ = − ∂ t A and B = ∇ × A using the Coulomb gauge, and neglecting longitudinal electric field in thefar-field region, we get (cid:104) : E µ ( r t ) B ν ( r (cid:48) t (cid:48) ) : (cid:105) t (cid:48) → t = ε νδξ (cid:16) − ∂∂x (cid:48) δ (cid:17) (cid:90) + ∞−∞ dω π (cid:126) ω (cid:104) D >µξ ( r , r (cid:48) , ω ) + D ( N ) ξµ ( r (cid:48) , r , − ω ) − D ( A ) µξ ( r , r (cid:48) , ω ) (cid:105) . (S5)Using the relations (see Ref.[S3]) D ( N ) µν ( r , r (cid:48) , ω ) = θ ( − ω ) D >µν ( r , r (cid:48) , ω ) , D ( A ) µν ( r , r (cid:48) , ω ) = θ ( ω ) D >µν ( r , r (cid:48) , ω ) , (S6)we get from Eq. (S5) (cid:104) : E µ ( r t ) B ν ( r (cid:48) t (cid:48) ) : (cid:105) t (cid:48) → t = ε νδξ (cid:16) − ∂∂x (cid:48) δ (cid:17)(cid:34) (cid:90) −∞ dω π (cid:126) ωD >µξ ( r , r (cid:48) , ω ) + (cid:90) + ∞ dω π (cid:126) ωD >ξµ ( r (cid:48) , r , − ω ) (cid:35) . (S7)Applying the relations D >µν ( r , r (cid:48) , ω ) = D <νµ ( r (cid:48) , r , − ω ) and D <µν ( r , r (cid:48) , ω ) = − (cid:2) D <νµ ( r (cid:48) , r , ω ) (cid:3) ∗ , we write Eq. (S7) as (cid:104) : E µ ( r t ) B ν ( r (cid:48) t (cid:48) ) : (cid:105) t (cid:48) → t = ε νδξ (cid:90) ∞ dω π (cid:126) ω (cid:20) − ∂∂x (cid:48) δ D <µξ ( r , r (cid:48) , ω ) (cid:21) . (S8)Bringing Eq. (S8) into Eq. (5), we get Eq. (7). Similarly, we can get Eqs. (8a) and (8b) in the main text.We write the free space GF in Eq. (S2) and its real space partial derivative to order O (1 /r ) as follows d rµν ( r, ω ) = C e qr r ( δ µν − ˆ R µ ˆ R ν ) − C e qr r q ( δ µν − R µ ˆ R ν ) + O (1 /r ) , (S9) ∂∂x γ d rµν ( r, ω ) = C e qr r q ˆ R γ ( δ µν − ˆ R µ ˆ R ν ) − C e qr r (cid:104) R γ ( δ µν − R µ ˆ R ν ) + δ γµ ˆ R ν + δ γν ˆ R µ (cid:105) + O (1 /r ) , (S10)with C = − / (4 πε c ) and ˆ R µ = x µ /r . Noting the full GF in Eq. (S1), we calculate the partial derivative of the fullGF in Eq. (7) as below ∂∂x (cid:48) δ D <νξ ( r , r (cid:48) , ω ) = d rνµ ( r , r i , ω ) χ <µ µ ( r i , r j , ω ) ∂∂x (cid:48) δ d aµ ξ ( r j , r (cid:48) , ω ) , (S11)where Einstein summation rule is assumed. Using the relation d a = ( d r ) † , and making the monopole approximation d rµν ( r , r i , ω ) = d rµν ( r − r i , ω ) ≈ d rµν ( r, ω ), we write Eq. (S11) as ∂∂x (cid:48) δ D <νξ ( r , r (cid:48) , ω ) ≈ d rνµ ( r, ω ) (cid:88) ij χ <µ µ ( r i , r j , ω ) ∂∂x (cid:48) δ (cid:2) d rξµ ( r (cid:48) , ω ) (cid:3) ∗ . (S12)Let’s first focus on the calculation of the energy radiation. We write the spherical surface integration of energyradiation by Eq. (10) as solid angle integration P = (cid:90) d Ω R S · ˆ R , (S13)where Ω is the solid angle of the spherical surface. In the large distance limit, R S · ˆ R in Eq. (S13) should be aquantity independent of the distance R due to conservation of the radiation energy. Thus, calculation of the Poyntingvector in Eq. (7) should be kept to order O (1 /r ), which implies that we only need to keep terms of order O (1 /r ) inEq. (S9) and Eq. (S10) for calculating the energy radiation. Eq. (S12) can thus be calculated as (cid:20) ∂∂x (cid:48) δ D <νξ ( r , r (cid:48) , ω ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) r = r (cid:48) ≈ C e qr r ( δ νµ − ˆ R ν ˆ R µ ) (cid:88) ij χ <µ µ ( r i , r j , ω ) C e − qr r ( − q ) ˆ R δ ( δ ξµ − ˆ R ξ ˆ R µ ) . (S14)Using Eq. (S14) to calculate Eq. (7) and bringing the result to Eq. (S13), we get the energy radiation power as P = − (cid:90) ∞ dω π (cid:126) ω π ε c (cid:90) d Ω( δ µν − ˆ R µ ˆ R ν ) (cid:88) ij Im (cid:2) χ <µν ( r i , r j , ω ) (cid:3) . (S15)We can make the approximation χ <µν ( r i , r j , ω ) ≈ Π <µν ( r i , r j , ω ) considering that the screening effect due to theradiation field is very small and can be neglected. Performing the solid angle integration in Eq. (S15), we get P = − (cid:90) ∞ dω π (cid:126) ω πε c Im (cid:2) Π tot ,<µµ ( ω ) (cid:3) , (S16)with Π tot ,<µν ( ω ) = (cid:80) ij Π <µν ( r i , r j , ω ). This is Eq. (12) in the main text.For the AM radiation, we write the surface integration in Eq. (11) as the solid angle integration d L dt = (cid:90) d Ω R (cid:104) : ˆ R × ←→ T : (cid:105) · ˆ R . (S17)In the large distance limit, R (cid:104) : ˆ R × ←→ T : (cid:105) · ˆ R in Eq. (S17) should be independent of R due to angular momentumconservation. This requires that the full GF and its partial derivative in Eq. (8a) and Eq. (8b) should be calculatedin the order of O (1 /r ). We need to keep the terms of order O (1 /r ) as well as O (1 /r ) in Eq. (S9) and Eq. (S10) tocalculate Eq. (8a) and Eq. (8b). The calculation is similar to that for the energy radiation. Bringing Eq. (8a) andEq. (8b) to Eq. (S17), we get the AM radiation rate as d L dt = (cid:90) ∞ dω π (cid:126) ω π ε c (cid:90) d ΩRe (cid:2) ˆ R × Π tot ,< ( ω ) · ˆ R (cid:3) . (S18)Performing the solid angle integration in Eq. (S18), we obtain Eq. (13) in the main text. THE SELECTION RULES FOR THE TIGHT-BINDING BENZENE MOLECULE In this section, we give some details of the derivation of the selection rules for the AM radiation of the TBbenzene molecule in the main text. To get an intuitive physical picture from the viewpoint that a photon is emittedvia a radiative transition by electrons from a state with higher energy to a state with lower energy, we make thetransformation from real space to mode space to analyze the radiation process. For the calculation of Eq. (12) andEq. (13) in the main text, the transformation to mode space is implemented asΠ tot ,<µν ( ω ) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) Tr (cid:104) M µ g < ( E ) M ν g > ( E − (cid:126) ω ) (cid:105) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) Tr (cid:104) U † M µ U U † g < ( E ) U U † M ν U U † g > ( E − (cid:126) ω ) U (cid:105) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) Tr (cid:104) (cid:102) M µ (cid:101) g < ( E ) (cid:102) M ν (cid:101) g > ( E − (cid:126) ω ) (cid:105) . (S19)Here M µ is the electron-photon coupling matrix summed over the photon site index k . For the simplest case thatonly two eigenstates with energy E m and E n of the TB benzene molecule are connected to the two leads respectively(see Fig. 2(b) in the main text). Eq. (S19) can be written asΠ tot ,<µν ( ω ) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) (cid:102) M µmn (cid:101) g 0. The first term in Eq. (S20) is zero, so we can writeΠ tot ,<µµ ( ω ) ≈ − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) (cid:102) M µnm (cid:102) M νmn (cid:101) g 6. Noting M µij = iev µij = i e (cid:126) t ij ( r i − r j ) µ , U jm = 1 / √ e i πjm/ , we get (cid:102) M x = eat (cid:126) i − i − i i − i − i i − i 00 0 0 2 i − ii i , (cid:102) M y = eat (cid:126) − − − − − − , (S24)( (cid:102) M x ) T · (cid:102) M x = (cid:16) eat (cid:126) (cid:17) , ( (cid:102) M x ) T · (cid:102) M y = (cid:16) eat (cid:126) (cid:17) i − i − i i − i i − i i 00 0 0 − i ii − i . (S25)Here, we have defined an element-wise matrix multiplication (cid:104) ( (cid:102) M µ ) T · (cid:102) M ν (cid:105) mn = (cid:102) M µnm (cid:102) M νmn . From Eq. (S24) andEq. (S25), we get the relations (cid:102) M xnm (cid:102) M xmn = (cid:102) M ynm (cid:102) M ymn , (cid:102) M xnm (cid:102) M ymn = i ∆ mn (cid:102) M xnm (cid:102) M xmn . (S26)Here, ∆ mn = sgn( m − n ) if | m − n | = 1, and ∆ = − ∆ = 1, otherwise ∆ mn = 0. Taking the weak coupling limit¯Γ → P = ω mn πε c (cid:16) (cid:102) M xnm (cid:102) M xmn + (cid:102) M ynm (cid:102) M ymn (cid:17) , (S27) dL z dt = ω mn πε c (cid:16) − i (cid:102) M xnm (cid:102) M ymn + i (cid:102) M ynm (cid:102) M xmn (cid:17) , (S28)with (cid:126) ω mn = E m − E n . Bringing Eq. (S26) into Eq. (S27) and Eq. (S28), we get dL z /dtP = ∆ mn ω mn . (S29)Since every emitted photon carries an energy (cid:126) ω mn , the number of photons emitted per unit time is dN/dt = P/ ( (cid:126) ω mn ). The angular momentum carried by an emitted photon is ∆ L = dL z /dtdN/dt = ∆ mn (cid:126) . DERIVATION OF THE RESONANT EFFECT In this section, we give a derivation of Eq. (14) in the main text for the resonant effect of the AM radiation. Theresonant effect is supposed to involve the inelastic transitions l = 2 → l = 1 and l = − → l = − (cid:126) ω ≈ t , supported by the FLSF with the chemical potential bias | µ L − µ R | = 2 t [Fig. 3(a)] and theAM radiation spectrum with the large peak at (cid:126) ω = 2 t [Fig. 3(b)]. The interacting self-energy including these twoprocesses is written as Π tot ,<µν ( ω ) = − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) (cid:102) M µ (cid:101) g < ( E ) (cid:102) M ν (cid:101) g > ( E − (cid:126) ω ) − i (cid:126) (cid:90) + ∞−∞ dE π (cid:126) (cid:102) M µ (cid:101) g < ( E ) (cid:102) M ν (cid:101) g > ( E − (cid:126) ω ) . (S30)Noting the (cid:102) M µ matrix in Eq. (S24), we getΠ tot , 1, we write the GFs as (cid:20) (cid:101) g r ( E ) (cid:101) g r ( E ) (cid:101) g r ( E ) (cid:101) g r ( E ) (cid:21) = (cid:34) E + t − (cid:2)(cid:101) Σ r leads (cid:3) − (cid:2)(cid:101) Σ r leads (cid:3) − (cid:2)(cid:101) Σ r leads (cid:3) E + t − (cid:2)(cid:101) Σ r leads (cid:3) (cid:35) − , (S32)and (cid:20) (cid:101) g < (cid:101) g < (cid:101) g < (cid:101) g < (cid:21) = (cid:20) (cid:101) g r (cid:101) g r (cid:101) g r (cid:101) g r (cid:21) (cid:34) (cid:2)(cid:101) Σ < leads (cid:3) (cid:2)(cid:101) Σ < leads (cid:3) (cid:2)(cid:101) Σ < leads (cid:3) (cid:2)(cid:101) Σ < leads (cid:3) (cid:35) (cid:20) (cid:101) g a (cid:101) g a (cid:101) g a (cid:101) g a (cid:21) . (S33)We consider the leads are coupled to the benzene molecule in the ortho position. The retarded lead self-energyis (cid:2) Σ r leads (cid:3) ij = − i Γ2 δ ij ( δ i + δ i ). The lesser and greater self-energies are (cid:2) Σ < leads (cid:3) = if L Γ, (cid:2) Σ < leads (cid:3) = if R Γ, (cid:2) Σ > leads (cid:3) = i ( − f L )Γ, (cid:2) Σ > leads (cid:3) = i ( − f R )Γ. The mode space self-energies can be obtained by (cid:101) Σ r,<,> leads = U † Σ r,<,> leads U . We get from Eq. (S32) (cid:101) g r ( E ) = 12 (cid:104) E + t + i Γ / 12 + 1 E + t + i Γ / (cid:105) ≈ E + t + i Γ / , (S34) (cid:101) g r ( E ) = 12 (cid:104) E + t + i Γ / − E + t + i Γ / (cid:105) ≈ i Γ / E + t + i Γ / , (S35)and (cid:101) g r ( E ) = (cid:101) g r ( E ). Similarly, we get (cid:101) g r ( E ) ≈ E − t + i Γ / , (S36) (cid:101) g r ( E ) ≈ i Γ / E − t + i Γ / , (S37)and (cid:101) g r ( E ) = (cid:101) g r ( E ). We obtain (cid:101) g < ( E ) (cid:101) g > ( E − (cid:126) ω ) − (cid:101) g < ( E ) (cid:101) g > ( E − (cid:126) ω ) ≈ Γ √ (cid:104)(cid:0)(cid:101) g r ( E ) (cid:1) ∗ (cid:101) g r ( E ) (cid:105)(cid:12)(cid:12)(cid:101) g r ( E − (cid:126) ω ) (cid:12)(cid:12) (cid:104) f L ( E ) − f R ( E ) (cid:105)(cid:104) − f L ( E − (cid:126) ω ) + f R ( E − (cid:126) ω ) (cid:105) − Γ √ (cid:104)(cid:0)(cid:101) g r ( E − (cid:126) ω ) (cid:1) ∗ (cid:101) g r ( E − (cid:126) ω ) (cid:105)(cid:12)(cid:12)(cid:101) g r ( E ) (cid:12)(cid:12) (cid:104) f L ( E − (cid:126) ω ) − f R ( E − (cid:126) ω ) (cid:105)(cid:104) f L ( E ) + f R ( E ) (cid:105) . (S38)In getting Eq. (S38), we have omitted some higher order terms considering that we may take the cross correlationsof degenerate states (cid:101) g r and (cid:101) g r as perturbations. Calculating Eq. (S31) using Eq. (S38) and bringing the result toEq. (13), we get the AM radiation dL z dt ≈ J (cid:34) θ ( − t − µ R ) (Γ / ( µ L − t ) + (Γ / − θ ( − t − µ L ) (Γ / ( µ R − t ) + (Γ / + θ ( µ R − t ) (Γ / ( µ L + t ) + (Γ / − θ ( µ L − t ) (Γ / ( µ R + t ) + (Γ / (cid:35) , (S39)with J = √ π tα ( v /c ) and α = e / (4 πε (cid:126) c ). In getting Eq. (S39), we have used zero temperature limit for the Fermifunctions of the two leads and assumed the lead coupling is weak, i.e., Γ (cid:28) t . Taking µ R = 4 eV > t in Eq. (S39), weget Eq. (14) in the main text.[S1] O. Keller, Quantum theory of near-field electrodynamics , (Springer, Berlin, Germany, 2012).[S2] K. Kaasbjerg and A. Nitzan, Phys. Rev. Lett. , 126803 (2015).[S3] M. Janowicz, D. Redding, and M. Holthaus, Phys. Rev. A , 043823 (2003).[S4] M. Katoh, M. Fujimoto, H. Kawaguchi, K. Tsuchiya, K. Ohmi, T. Kaneyasu, Y. Taira, M. Hosaka, A. Mochihashi,and Y. Takashima, Phys. Rev. Lett.118