Ultrafast microscopy of a twisted plasmonic spin skyrmion

Topology engages ideas within solid state materials spanning the nanometer to cosmic scales from subatomic particles to cosmology. Topology can be intrinsic to the lattice and consequent electronic structure, or be imposed on them by intense electromagnetic fields. Skyrmions are known as stable magnetic quasiparticles of interest to fundamental physics, and with potential applications to magnetic data storage and quantum computing. Research on the photonic and plasmonic skyrmions, however, is at its infancy. We present the first ultrafast microscopy imaging study of a plasmonic spin skyrmion generated by the spin–orbit interaction of light. We reveal by experiment and theory the nanofemto plasmonic vortex fields and their topological skyrmion spin textures that can dress topologically trivial materials, drive magnetoelectric phenomena, and probe axion physics.

Topological solitons are singular field distributions broadly relevant to the physics and chemistry of structure formation, charge-parity violation, phase transitions relevant to cosmological evolution, nucleons, and condensed matter.^1–7 Skyrmion solitons were first proposed as homotopically nontrivial differential vector field theoretic descriptors of the structure of nucleons⁸ and later were discovered in condensed mater contexts as more accessible homotopes, such as stable field and spin textures in chiral magnets,^9,10 Bose–Einstein condensates,^11,12 liquid crystals,¹³ free space vectorial optical fields,^14–17 propagating acoustic waves,¹⁰⁹ and plasmon^7,18,19 as well as phonon polaritonic fields.²⁰ 

The discovered plasmonic field skyrmion textures form with locally linearly polarized surface plasmon polariton (SPP) standing wave electric fields oscillating at PHz frequency,¹⁸ such that their topological textures average to zero within each optical cycle and they do not possess a chiral character.¹⁹ Photonic spin angular momentum (SAM), which defines the chiral quantum interactions of light with matter, can also be engineered to possess skyrmion textures,^7,21 whose topological properties, however, are expected to be stable over the light pulse duration, as has been demonstrated in the half-skyrmion (meron) textures.^22,23 Such topological spin textures, which have escaped scrutiny until methods that can probe them on a deep subwavelength scale had been developed, are formed within plasmonic vortices carrying orbital angular momentum (OAM), despite that structured OAM light has been known for more than 30 years.⁶ 

Unlike ground state skyrmions, such as those in magnetic materials that form through spin-spin interactions on picosecond time scales,²⁴ plasmonic spin skyrmions exist and can be investigated in ultrafast dynamical contexts where they can impress their properties on matter. When they are created and persist on few femtosecond, or longer timescales of the generating pulse duration, such topological fields can perform Poincaré engineering of matter by impressing their topological field and spin textures on proximate materials in their near field;^25,26 for example, they can induce topological phase transitions,^10,22,24 nonreciprocal magnetoelectric responses,^27,28 and valley selective interactions in 2D semiconductors.^29–32 They can also control topological quasiparticle functions in trivial materials^26,33,34 and expand repertoires of optoelectronics, single molecule spectroscopy,^35–37 and topological photonics.^38–40 

Quantitative experimental information on the plasmonic spin skyrmion nanoscale textures, their skyrmion boundaries, and integer topological charges, however, is yet to be obtained because of the challenge in the imaging of the nanofemto SPP fields. Moreover, demonstration of their dynamic stability on femtosecond scales, as well as the associated intense nanoscale electromagnetic fields, which are of particular interest to chiral, nonreciprocal light–mater interactions, are lacking. Here, we directly address these questions by reporting an analytical method for the generation of twisted skyrmion spin textures and their associated chiral electromagnetic fields at the core of plasmonic vortices, and their quantitative visualization through coherent nanofemto imaging by interferometric time-resolved photoemission electron microscopy (ITR-PEEM).^41,42 Our analytical model finds that the photonic spin–orbit interaction (SOI) within plasmonic vortices with quantized OAM of $l$ = ±1 molds the plasmonic SAM into a chiral skyrmion texture in the vacuum side, and a skyrmion counterpart of opposite topological charge in the metal side, whose cores host co-gyrating, anti-parallel in-plane E and B fields (E·B focus) that can drive magnetoelectric phenomena, such as chiral axion electrodynamics.^43–45 Together, the skyrmion pair forms a hedgehog spin texture within the primary vortex. The boundaries of the plasmonic spin skyrmion are defined by contours of linear polarization known as L-line singularities,^46,47 where atomically sharp superposition cat states of the plasmonic up- and down-spins are observed, which invite applications in quantum information processing at the atomic scale.^48,49 By recording the local field induced nonlinear photoelectrons using ITR-PEEM with nanometer resolution and attosecond scanning accuracy,^42,50–52 we image the dynamical vectorial vortex fields, based on which we extract the plasmonic spin skyrmion textures,⁵³ skyrmion boundaries, and topological charges, which are in quantitative agreement with our analytical model. The skyrmion topological character is also demonstrated to be stable longer than ∼10 optical cycles; thus, it can potentially function as information bits,^54–56 aid in Fourier optical processing,⁵⁷ transiently break the time reversal symmetry (TRS) to actuate coherent magnetoelectric devices,^44,58–61 and act as active nodes in quantum information storage and processing.⁶² 

SPPs are bound evanescent fields that propagate at metal/dielectric interfaces.⁶³ When circularly polarized light with helicity σ illuminates a sharp circular slit structure (see the supplementary material Fig. S1), in an otherwise flat polycrystalline Ag film,²² an SPP vortex with OAM of $l=σ$ is generated.⁶⁴ Equations (1a)–(1c) describe the vectorial field components (⁠$Ez$⁠, and $E∥=Er+Eθ$⁠) of a plasmonic vortex carrying OAM of $l$ that forms a phase singularity in $Ez$ at the vortex core,^64–66 where $r$ and $θ$ represent the chosen cylindrical coordinates (see the supplementary material),

In Eq. (1), $0+$ refers to the dielectric side, $E0$ is the field amplitude which is normalized to 1 afterward, $ε0 (εd)$ is the absolute (relative) bulk permittivity of the vacuum, $i=−1$⁠, $ωSPP$ is the SPP frequency, $kSPP$ and $κ$ are the SPP wave vectors in-plane and normal to the interface, $R$ is the radius of the coupling structure, and $Jv$ is a Bessel function of order $v$ that is defined by OAM of the orbiting $Ez$ component, i.e., $v$ = $l$ (see the supplementary material, and Fig. S1). Note that an imaginary part appears in Eq. (1b), which accounts for the finite $R$⁠, and causes the in-plane spin to twist, as will be shown later. In the following, we only consider the z = 0⁺ SPP field components at the interface, because they generate the two-photon photoemission (2PP) signal that we image.

Pseudo color plots of the distribution of the field components from Eqs. (1a)–(1c) are shown in Figs. 1(a)–1(c), for an $l$= −1 vortex with $R=2λSPP$⁠. The $|Ez|$ component has maximum at $r1≈0.29λSPP≈150$ nm for a plasmon wavelength $λSPP$=530 nm, to which we refer hereafter as the primary vortex ring [small dashed circle in Fig. 1(a)]. The superposed arrows in Fig. 1(a) indicate the in-plane Poynting vectors, corresponding to a counterclockwise (CCW) SPP propagation; the time evolution of the corresponding vortex field components are shown in Fig. S2. Moreover, while $Ez$ vanishes at the vortex core where the phase is singular, a feature known as superoscillation,⁶⁷ the $|Er|$ and $|Eθ|$ have a strong focus, with $|Er|$ passing through zero at $r1$⁠, and $|Eθ|$ diminishing more slowly with $r$⁠. The minimum of $|Eθ|$ occurs at $r2≈0.61λSPP$⁠.

A novel aspect of vortex fields that emerges from this analytical model, and which is likely to find important applications in magnetoelectric interactions and topology, is that vortex cores are the focus of E$∥$B fields (Fig. S2). Although we do not experimentally image the plasmonic B field of the SPPs because of the lack of dipole interaction that can generate 2PP signal, their in-plane (⁠$H∥=Hr+Hθ$⁠) components also form a focus at the vortex core, which resembles $E∥$⁠, but have the opposite sign [SI Eq. (S4)]. Therefore, the surface parallel $E∥$ and $B∥=μH∥$ fields combine as a co-gyrating antiparallel pair in a tight focus at the vortex core. Like the photonic spin texture, the co-gyrating antiparallel relationship persists as long as the vortex is established. Due to the strong focusing of the in-plane fields, which drives the nonlinear photoemission, we expect vortex cores to be a focus of intense magnetoelectric coupling defined by the pseudoscalar $E∥·B∥$ product. The calculated distribution of $E∥·B∥$ from the analytical model is shown in Fig. 1(d), where the anti-parallel $E∥$ and $B∥$ fields are focused on the <$λSPP$/2 length scale; the nanofemto distribution of the counterclockwise co-gyrating in-plane field vectors is shown in the SI movie 1 (Multimedia view). Note that the $E∥$ and $B∥$ fields become parallel and co-gyrate in a clockwise direction when the sign of vortex OAM is reversed.

The (anti)parallel context of the $E∥$ and $B∥$ fields at the vortex core is a consequence of circularly polarized light interacting by the circular coupling structure; it is thus a profoundly simple construct for interrogating magnetoelectric phenomena by ITR-PEEM and other near-field microscopic techniques. Moreover, such (anti)parallel E and B fields are a property of vortices composed of surface evanescent waves, in general, where the associated transverse field forms a singularity at regions of the strong E $∥$ B focus due to the cycloidal field polarization in the meridional plane. Therefore, such fields can also form in 2D semimetals/semiconductors where plasmon frequency is in the THz regime, or in 2D insulators in the form of phonon polaritons, etc.^20,68 Consequently, the E $∥$ B focus can potentially drive magnetoelectric effect in a broad frequency regime, such as magnetoelectric coupling in metallic plasmonic nanoparticle dimers, where it drives coherently the light-nanostructure orientation-dependent magnetic and electric dipole chiral Raman scattering,⁶⁹ or excitation of magnetic moments parallel to the driving electric fields in bianisotropic metamaterials.^38,70–72 Finally, the strong E $∥$ B may assist in the study of dynamical axion field $θ$ through the topological magnetoelectric coupling $θE·B,$^43,73–76 because it locally enhances the electromagnetic field strengths, breaks the TRS on subwavelength scale, and its frequency response can be tuned continuously to fall within the bandgap of magnetically doped topological insulators.

Having defined the vectorial SPP field vortex components in Eq. (1) and Fig. 1 for $l$= −1, we use them to calculate the associated SAM (S) profile, defined as^77,78

where $ω$ is frequency of the time-dependent electric $E$ and magnetic $H$ fields, and the asterisks indicate complex conjugation. According to the derived $E$ and $H$ field components (supplementary material), the vacuum SAM components for $l$ = −1 are

The $Sz$ component is discontinuous, and on the metal side (⁠$0−)$ has an equivalent form, but the $εd→εm$ (metallic permittivity) substitution being below the plasma frequency, reverses its sign.

Based on Eq. (3), we plot the normalized vacuum-side skyrmion SAM texture, $Sn=S|S|$ in Fig. 2(a). In this texture, vacuum SAM points up (+z) at the vortex core, aligning it with the SAM direction of the excitation light; moving away from the core, it continuously winds to the surface plane pointing away from the core at $r1$⁠, and finally down at its periphery at $r2$⁠,^79,80 as is illustrated in Fig. 2(b). At $r2$⁠, $Sz,n$⁠, the z-component $Sn$⁠, has a π-phase reversal [Figs. 2(b) and 3(b)], such that the $Sz,n$ vectors on each side of $r2$ are surface normal, but point in the opposite directions. We note, however, that the unnormalized spin, $Sz$⁠, is continuous as it smoothly reverses direction upon crossing through zero at $r2$⁠.⁷ The anti-parallel $Sz,n$ vectors across $r2$ reflect the counterrotation of the in-plane fields that originate from the evanescent properties of SPPs. The superposition of the oppositely circulating in-plane fields at $r2$ necessitates a transition through linear polarization, which is a chirality discontinuity catastrophe that occurs on the screening length scale where chiral interactions can achieve deep subwavelength contrast,^81–83 and the electron transport direction is reversed.

The in-plane component of the normalized SAM, $S∥,n=Sr,n+Sθ,n$ [arrows in Fig. 2(c)] has a phase that is spatially shifted from $Sz$ by ∼$14λSPP$⁠. Consequently, its amplitude is zero at the vortex core and at $r2$⁠, but reaches maximum at $r1$⁠, as seen by the pseudoscalar plot of $|S∥,n|$ in Fig. 2(c). Note that there is also a finite $Sθ,n$ proportional to $R$⁠, which is evident as a twist of $S∥,n$ from the radial direction in Fig. 2(c). It appears because $r$ is not negligible with respect to $R$ in deriving the Bessel function orders. Therefore, the finite $Sθ$ turns the SAM from a Néel- toward a twisted Bloch-type skyrmion texture when $r∼R$⁠.⁸⁰ The $R$-dependent SAM twist angle is quantified in supplementary material. Interestingly, plasmonic spin skyrmion textures possessing similar twisting characters have also been theoretically demonstrated to occur in chiral materials in between metallic films, but originating from a distinct mechanism.⁸⁴ 

Similar to the interference of spin waves in magnetic skyrmion lattices, where a skyrmion unit is defined by a soliton solution of interfering waves with a defined wavelength,^85,86 we define a radius dependent skyrmion boundary based on the SPP interferences governed by the zeros of Bessel functions, and the Poincaré sphere of light polarization [Fig. 2(d)].^87,88 The first boundary is marked by the dashed circle of radius $r2$ in Figs. 2(b) and 2(e). Such boundary is intrinsically associated with the second L-line polarization singularity fringe of the in-plane SPP fields. The L-line singularity⁴⁶ defines the interface where the $E∥$ field is linearly polarized, and consequently, the out-of-plane SAM passes through $Sz=0$⁠. Mapping the SPP polarization onto a Poincaré polarization sphere places the L-line at the equator, along which the linear $E∥$ field polarization rotates by 2π [Fig. 2(d)]. From the Poincaré sphere perspective, the surface-normal positive and negative $Sz$ components at the vortex core and at $r2$ are the north and south poles, respectively, representing circular polarization, or C-point singularities;^46,47 the upper hemisphere bounded by the equator (⁠$r1$⁠) encloses a half-skyrmion, or a meron [Fig. 2(d)].^22,23,89 Note that at L-lines, the out-of-plane SAM vanishes, while the in-plane SAM is maximum due to the momentum-SAM locking. We numerically compute the L-line and C-point distributions by calculating the ellipticity of the $E∥$ fields, as shown in Fig. 2(e). The minimum ellipticity of 1 defines the C-point singularity at the vortex core [the center in Fig. 2(e)].^46,47 Away from the C-point, the SPP field becomes more elliptical, passing through linear at the first L-line singularity at $r1$ [$Sz→0$ in Fig. 2(b)].

To evaluate the topological charge of the plasmonic SAM texture, we calculate the charge density, $D$⁠, as defined by⁹⁰ 

where the structure factor, $Sn$⁠, is the normalized SAM at each pixel in the Cartesian coordinates. The calculated areal distribution of $D$ for the $l$ = −1 vortex is shown in Fig. 2(f). Integration of D within the area defined by the second L-line [Fig. 2(e)] gives an integer topological charge of N = +1 for the $l$ = −1 plasmonic vortex, confirming its skyrmion topological character.²⁴ For an $l$ = +1 SPP vortex excited with right circularly polarized light, the topological charge is N = −1, due to the change in sign of the SAM pseudovectors.

The plasmonic skyrmion spin texture, defined here by the L-line and C-point analysis, establishes its homotopy to that of a twisted magnetic skyrmion,^46,47 and defines a precise method to locate both their dominant density (⁠$r1$⁠) and boundary (⁠$r2$⁠). Finally, because the plasmonic spin winding recurs along the radial direction as described by Bessel functions, the analysis of the topological texture based on Eq. (4) can be extended outward from the vortex core.

A notable aspect of plasmonic skyrmion SAM textures is that they form on both sides of the metal/vacuum interface, but with the opposite topological number according to the photonic spin topology across the interface described by a strong topological insulator ℤ₄ index;^91–93 the interface-normal cross-section of the normalized SAM vectors $Sn$ is shown in Fig. 3(a), and the vectorial spin components on vacuum and metal side of the interface are shown in Figs. 3(b) and 3(c), respectively. The $Sn$ rotates by 2π with respect to the vortex core within $r1$⁠, because ε, and thus $Sz∼εE*×Ez$⁠, changes sign on an atomically sharp scale across the interface;⁹¹ the in-plane SAM, which is locked to the CCW circulating k_SPP-vector, however, remains parallel across the interface. Within $r1$⁠, the combined spin texture of the two half-skyrmions form a hedgehog texture, which is homotopic to Berry curvature monopoles in topological materials,^74–76,94 magnetic monopole fields for electron vortex beam generation,^95,96 and spin monopoles in Bose–Einstein condensates⁹⁷ and photonic systems.⁹⁸ We note, however, the rate of spin texture winding is asymmetric across the interface, which is due to the joint electric and magnetic contributions to the total SAM (see supplementary material).

Unlike a macroscopic TRS breaking induced by a magnetic phase transition, for plasmonic vortices, the TRS is broken on λ_SPP length and the optical pulse time scales by a gyrating electron current that induces an inverse Faraday effect,⁶¹ and the concomitant surface normal plasmonic SAM. Therefore, plasmonic spin hedgehog can potentially induce spin-Hall effect in transport of electrons on the ultrafast timescale.^99,100 Moreover, it may also act as sources of photoelectrons carrying OAM from the SPP field that is of interest for electron microscopy of magnetic materials.^96,101 We note that, compared to the in-plane skyrmion-anti-skyrmion pairs in magnetic materials,^102,103 the plasmonic skyrmion spin textures are more tightly confined at the quasi-2D metal/dielectric interfaces defined by the plasmon field decay.

To validate the SAM skyrmion texture within an |$l$| = 1 vortex core and its temporal structure, we record movies of the vectorial SPP vortex field dynamics by imaging the spatial distribution and temporal evolution of the 2PP signal by the ITR-PEEM method.^42,52 Figure 4(a) shows a schematic view of the ITR-PEEM experiment where two circularly polarized pulses of mutual delay τ illuminate a polycrystalline Ag film with a fabricated ring SPP coupling structure. The coupling structure causes the SPP field focusing with an orbital angular momentum defined both by the geometric charge of the coupling structure and spin of the generating optical field. The as-recorded, one pulse, photoemission image obtained by left circularly polarized light excitation is shown in Fig. 4(b); the 2PP signal comes to an intense focus within a circular region of ∼$λSPP$ diameter centered on the vortex core. Excitation by light of the opposite handedness generates superficially the same 2PP distribution (see supplementary material). All SPP field components contribute to the total 2PP signal near the vortex core; this superposition of the in- and out-of-plane field-induced signals shrouds the SAM skyrmion, because one cannot immediately separate them.

The separate SPP field component-induced 2PP signals are obtained in an ITR-PEEM experiment by Fourier filtering the time-domain ITR-PEEM image sequence (see supplementary material).²⁵ This generates a phase dependent movie of the $|Er|$ SPP field component oscillating at the laser driving frequency 1ω_L (SI movie 2). (Multimedia view) Fig. 4(c) reports the Fourier filtered image of the $|Er|$ component at τ ∼ 23.1 fs (τ = 0 fs is defined by the pump-probe pulses' temporal overlap), where the image symmetrization and smoothing is applied as described in supplementary material. We further derive the $|Ez|$ and $|Eθ|$ components from the $Er$ distribution,⁵³ by applying the Maxwell's equations and plot them in Figs. 4(d) and 4(e). The dashed circles in Figs. 4(c)–4(e) mark the location of the vortex ring at $r1≈150$ nm, confirming agreement between the experiment and the analytical model [Eq. (2)]. Furthermore, from the experimental fields shown in Figs. 4(c)–4(e), we construct the vectorial spin textures based on the corresponding Bessel function distributions. In Fig. 5(a), we show the derived distributions of $Sz,n$⁠, and in Fig. 5(b) the magnitude and the directions of the normalized in-plane SAM components $S||,n$ by a pseudo color map and arrows, respectively, which are in excellent agreement with the predicted distributions in Figs. 2(b) and 2(c). The observed continuous winding of $Sz,n$ from pointing up at the center to down at ∼$r2$⁠, and in particular, with the distinct twist of $S||,n$ unambiguously confirm the generation of twisted skyrmion plasmonic spin texture at the $l$ = −1 vortex. To avoid confusion, however, we do not directly image the SPP SAM, but our experiment enables to extract the vectorial fields, and based on the Maxwell's equations, to transform them into their consequent spin distributions.²⁵ 

Moreover, the spin, and specifically, the inverse of the $|Sz|$⁠, defines the in-plane SPP field ellipticity (L-line) distribution at the plasmonic vortex. Figure 5(c) shows the experimental ellipticity distribution of the $l$= −1 vortex, where distinct sharp contrast appears at both at $r1$ and $r2$⁠, because the ellipticity becomes infinite and the SAM points in-plane. Such analysis locates the L-lines from the derived SAM textures. This L-line distribution agrees remarkably well with the calculated one in Fig. 2(e), to experimentally locate the stable skyrmion texture that is bounded by the second L-line from the vortex core. Conclusive experimental evidence of the topological texture assignment is shown in Fig. 5(d), where the L-line distribution (blue) is plotted, along with the integrated topological charge density (red) as a function the integration radius. Horizontal black dashed lines indicating the location of integer topological charges are also shown intersecting with the even-numbered L-lines, in agreement with the numerical calculations shown by corresponding dashed curves in Fig. 5(d). Therefore, the even-numbered L-lines naturally define boundaries separating quantized integer increments in the topological charge that confirm the homotopy of the SAM skyrmion textures. The odd L-lines, however, mark the half-skyrmion, or meron textures, where the quasiparticle density [Fig. 2(f)], and thus, the topological charge varies the most rapidly.

Finally, with the capability of resolving SPP vortex field dynamics supporting the spin skyrmion quasiparticle with nanofemto precision, we demonstrate that such quasiparticle is stable over the SPP pulse duration. Figure 5(e) plots the spin skyrmion topological charge as a function of inter-pulse delay τ in unit of optical cycles; the integer topological charge is observed to stabilize after ∼20 fs when the SPP vortex is fully established, and then remains constant for over ∼10 optical cycles (∼18 fs) before it dissolves. Therefore, such dynamic spin quasiparticle, which breaks the TRS on deep subwavelength scale, can serve pulses of transient magnetization at THz frequencies, which can actuate the inverse Faraday effect or ultrafast spin-dependent electron transport switching. The temporal stability of topological charge is defined by the generation pulse width and thus can be set by illumination with circularly polarized fields of any duration. To exist in the absence of the generating electromagnetic field, as does the magnetic skyrmion, would require the magnetic monopole spin texture to be prepared by other means.

Our study unambiguously demonstrates, with nanofemto resolution, an ultrafast twisted plasmonic SAM skyrmion texture, its sub-diffraction limited polarization L-line singularity boundaries, and integer topological character, which are dynamically stable on ∼20 fs timescale of the generation pulse duration. In addition, our model also finds that (1) there exists an intense focus of the $E||·B||$ term that is of potential interest for nonlinear magnetoelectric couplings and axion field dynamics in quantum materials,^43,73,104 and (2) a pair of skyrmions of opposite charge are generated across the SPP supporting interface, where within the primary vortex ring the plasmonic spin can be considered as a monopole-like hedgehog texture, and can induce a transient spin-Hall effect.^99,100 The topological spin textures generated by the interaction of optical fields with OAM structures that induce polaritonic spin–orbit interaction is a general phenomenology of focusing of polaritonic fields that carry OAM.²⁰ Finally, because the reported plasmonic skyrmion is delivered by intense plasmonic fields concentrated on the λ_SPP scale, it is expected to imprint its topological spin texture and OAM on topologically trivial or nontrivial materials^105,106 on femtosecond timescale. Thus, the SAM skyrmions are ideal for provoking topological quasi-energy states by Floquet engineering, and texturing of spins by Poincaré engineering to achieve TRS breaking, and nonreciprocal physics in the condensed quantum matter.^{38,39,43,107,108}

See supplementary material for detailed derivations of the twisted plasmonic skyrmion textures, as well as the processing and analysis of the ITR-PEEM data. Figures S1–S6, Movies S1 and S2. (Multimedia view).

This research was supported by the NSF Center for Chemical Innovation on Chemistry at the Space-Time Limit Grant No. CHE-1414466 and ONR MURI Grant No. N00014-20-S-F003. Y.D. thanks R. Mong for consultation on the analytical theory; H.P. thanks V. Liu and V. A. Apkarian for discussion of magnetoelectric effects and axion physics.

The authors have no competing interest.

Y.D. performed the analytical modeling and experiments; Z.Z., A.G., and K.K. contributed to the nanofabrication, electromagnetic simulation, and laser operation; M.D. supervised the establishment of the ITR-PEEM experiment; A.K. advised on the ITR-PEEM measurements and data analysis; C.B.H. introduced and guided the research on topological plasmonics; H.P. supervised the research; Y.D. and H.P. wrote the manuscript.

The theoretical simulation results and experimental ITR-PEEM data are available from the authors upon request. The analytical model of topological skyrmion field textures is presented in the paper and the supplementary material.