Skyrmions are topologically protected quasiparticles characterized by three-dimensional (3D) topologically nontrivial spin textures into a two-dimensional (2D) plane via stereographic projection. The concept was originally proposed by nuclear physicist Tony Skyrme in the 1960s as soliton solutions in nonlinear field theory [1]. Since then, skyrmions and other quasiparticles have been widely predicted in various physical systems, including Bose–Einstein condensates [2], liquid crystals [3], superfluid [4], magnetic materials [5], and twistronics [6]. Particularly, magnetic skyrmions, owing to their topological stability and nanoscale feature size, are considered promising candidates for next-generation data storage [5]. Recently, skyrmionic structures have been introduced into optics [7–11], where topological skyrmionic textures are constructed by various degrees-of-freedom of electromagnetic (EM) fields, including the spin angular momentum (SAM) [8,10,12–23] and the electric-field vector [7,9,24–26] in confined fields, the Poynting vectors [27,28] and the Stokes vectors [29–40] in propagating fields, and the pseudo-spin vectors in nonlinear photonic crystals [41]. Similar to topological quasiparticles in magnetic materials [42–44], optical skyrmions can form either as isolated quasiparticles [11,28–39] or in lattice configurations [7,9–12,14–18]. For instance, in surface plasmon polaritons (SPPs) with continuous rotational symmetry, the spin-orbit coupling (SOC) of light leads to the formation of skyrmion-like spin quasiparticles [8], whereas in discrete rotational symmetries, spin skyrmion or meron lattices emerge, with their topologies constrained by both the rotational symmetry and the conservative property of total angular momentum (TAM) [13,20,21,45–48]. Optical quasiparticles show great potential for applications in optical communications [49,50] and topological quantum communications [51], and the concept has also been extended to other physical systems, such as water waves [52,53] and acoustic waves [54].

Previously, most research on topological quasiparticles in the optical near field has primarily focused on the construction of various novel types of optical topological textures at the metal/isotropic media interface, and the topological stability of the resulting textures remains unexplored. In particular, the evolution of optical topological lattices in anisotropic media has yet to be investigated. The anisotropy related to the lattice symmetry of crystals or artificial meta-atoms is widely present in modern nanophotonics, such as liquid crystal devices for topological storage and logic devices [55], and metamaterials and metasurfaces for chiral bound states in the continuum [56,57], and thus the investigation of the evolution of optical topological lattices in anisotropic media is critical for the future exploitation of applications.

In this study, we systematically investigated the formation, topological stability, and phase transitions of optical meron lattices at the metal/uniaxial crystal interface. We derived the explicit expressions for the electric and magnetic fields, as well as the SAM, of surface plane waves at the metal/uniaxial crystal interface. By arranging these surface plane waves in C4 rotational symmetry and considering the conservation of TAM, we constructed electric-field meron lattices and spin meron lattices, respectively. More importantly, our theoretical results demonstrate that varying the optical axis (OA) orientation of the uniaxial crystal can induce phase transitions in the electric-field meron lattice, from a topologically nontrivial to a topologically trivial state. In contrast, the skyrmion number of the spin meron lattice remains stable against such rotations, demonstrating its topological stability. This difference arises because in the instance of the electric-field meron lattices, the surface waves are no longer pure transverse magnetic (TM) SPP modes at the metal/uniaxial crystal interface, whereas for the spin meron lattices, the decomposition of TAM into orbital and spin components in an anisotropic medium is similar to that in an isotropic medium. The results provide a novel mechanism to achieve the topological stability and phase transitions of topological quasiparticles in optical systems and open new avenues for advanced spin-optics and topological photonics applications.

The optical system as illustrated in Fig. 1(a) consists of a semi-infinite uniaxial crystal occupying the region z>0 and an isotropic metal in the region z<0. The x and y axes lie in the interface between the two media (z=0), while the z axis is perpendicular to the interface. Silver is selected as the isotropic metal material with its dielectric constant set to εm=−18.294+0.481i at a wavelength of λ=633nm [58]; the imaginary part is neglected for simplicity. Without loss of generality, we first assume that the principal axes of the uniaxial crystal are aligned along the x, y, and z directions of the Cartesian coordinate system. The principal dielectric constants and refractive indices of the crystal are denoted as εx,εy,εz and nx,ny,nz, respectively. In this case, the permittivity tensor of the uniaxial crystal can be expressed as a diagonal matrix comprising the principal elements: diag[εx,εy,εz]=diag[ε0nx2,ε0ny2,ε0nz2], where ε0 is the permittivity in vacuum. The anisotropic medium is chosen as the negative uniaxial crystal material calcite (CaCO3) [59], with the ordinary refractive indices ny=nz=no=1.6557 and the extraordinary refractive index nx=ne=1.4849.

To derive the expressions of electric and magnetic fields of the hybrid surface waves and the propagation constants, we consider the case where the hybrid surface plane waves propagate along the x direction and rotate the OA orientation angle Φ, which is defined by the intersection angle between the OA and the x axis. In this instance, the corresponding permittivity tensor can be obtained by the rotational operator as (see Appendix A for details) ε↔=(ηxxηxy0ηyxηyy000ηzz),(1)where ηxx=ε0(nx2cos2Φ+ny2sin2Φ), ηyy=ε0(nx2sin2Φ+ny2cos2Φ), and ηxy=ηyx=ε0(nx2−ny2)sin2Φ/2. It can be observed that, when the OA is rotated, the horizontal components and the nondiagonal terms determined by the no,ne, and Φ of the permittivity tensor are changed correspondingly, which results in the variations of the propagation constants.

The metal/uniaxial crystal waveguide can support the hybrid SPP modes that include both TM and transverse electric (TE) polarization components [60,61]. As indicated above, we assume that the surface plane wave propagates along the x axis and decays exponentially in the z direction: Ez=E0exp(iβx−kzz) and Hz=H0exp(iβx−kzz), where β and ikz are the transverse (propagation constant) and longitudinal (attenuation coefficient) wave vector components. Based on Eq. (B7) in Appendix B, we obtain that the propagation constant in the hybrid surface mode is determined by β=(βnum+βnum2+4ηzzβden/ηxx)/2,(2)where there are βnum=βyy2+(ηzzβxx2−ω2ηxyηyxμ)/ηxx, βden=βxy2βyx2−βxx2βyy2 with βxx2=ω2ηxxμ+kz2,βyy2=ω2ηyyμ+kz2 and βxy2=βyx2=ω2ηxyμ=ω2ηyxμ. Here, ω represents the angular frequency of the EM field, and μ denotes the permeability. By applying the boundary conditions, the attenuation coefficient kz and the propagation constant β can be obtained. As shown in Fig. 1(b), the normalized parameter β/k0 reaches its maximum at Φ=45° and Φ=135°, and attains its minimum at Φ=90°, where k0 denotes the wavenumber in vacuum. Subsequently, the distributions of the electric and magnetic fields via the OA orientation Φ can be derived (see Appendix B for details).

Optical topological lattices can be constructed by the superimposition of multiple surface plane waves. However, it is worth noting that the fundamental mechanisms for generating the optical lattices are the rotational symmetry and the conservative property of TAM in optical systems [13,20,21]. By superimposing multiple surface plane waves, it is possible to approach the solutions of Maxwell’s equations under discrete rotational symmetries. Therein, optical topological lattices with strictly rotational and translational symmetries can be formed in C3,C4, and C6 rotational symmetries, whereas other polygonal symmetries, such as five-fold rotational symmetry, give rise to optical quasi-crystals, moiré lattices, and mixing skyrmions and merons topologies [45–48]. In systems with continuous rotational symmetry, the isolated skyrmions can be formed [8]. In C4 rotational symmetry, the out-of-plane electric-field component Ez of the hybrid surface wave can be derived by the Hertz potential, expressed as [62] Ez=Ψ4=E0∑n=14eilθneiβn(xcosθn+ysinθn)e−kzz,(3)where E0 is a constant, θn=nπ/2, and l represents the quantum number of TAM.

In the real optical systems for generating optical topological lattices, it is more convenient to fix the uniaxial crystal and rotate the excited surface waves. Although Eq. (2) considers the solution of the propagation constant when the surface wave propagates along the x axis direction and the OA of the uniaxial crystal is rotated, it can be transformed to the solution of the real optical systems by a rotation matrix along the z axis flexibly. Consequently, Ez can be reformulated as Ez=E0∑n=14eilϑneiβn(xcosϑn+ysinϑn)e−kzz,(4)where ϑn=nπ/2+Φ, and Φ plays a key role in tuning the optical topological texture. It can be observed from Eq. (4) that even in the presence of an anisotropic uniaxial crystal, the distributions of EM fields still exhibit periodicity with respect to the quantum number of TAM.

The topological property of the skyrmionic texture can be characterized by the skyrmion number, which is defined as follows [5]: Q=14π∬n·(∂n∂x×∂n∂y)dxdy.(5)

Here, n denotes the normalized unit vector, and the skyrmion number quantifies how many times the vectors wrap around a unit sphere. Specifically, Q=±1 corresponds to skyrmions, in which the entire unit sphere is covered once, whereas Q=±1/2 corresponds to merons, in which only a hemisphere is covered [43,44].

In this section, we first analyze the case l=0 and the propagating directions of surface waves align with the principal axes (i.e., Φ=0°). In this instance, the four surface plane waves form the standing wave fields along the x and y directions. The standing wave fields interfere at the metal/uniaxial crystal interface, thereby generating the electric-field meron topological lattice. Figure 2(a) shows the real part of the out-of-plane electric-field component Re(Ez), which exhibits a C4 symmetry. The corresponding 3D electric-field vector structure is shown in Fig. 2(b). Within each unit cell, the electric-field vector varies from a central “up” or “down” state to an edge flat state, which manifests as a Néel-type meron texture with alternating polarity. According to Eq. (5) and n=Re(E)/|E|, we use the breadth-first search (BFS) algorithm [63] to identify the boundaries where nz=0, and calculated the skyrmion number within the unit cell marked by a black square to be Q=0.493. The skyrmion number can approach 0.5 by increasing the number of sample points further.

To evaluate the stability of the electric-field meron lattices with respect to the OA orientation of the uniaxial crystal, we rotate the four surface plane waves and compute the skyrmion number of the central unit cell under varying Φ, as illustrated in Fig. 2(c). The skyrmion number is found to oscillate as the Φ varies from 0° to 180°, and it reaches 0.5 only at discrete angles: Φ=0°, 45°, 90°, 135°, and 180°. The results reveal that the topological properties of the optical lattices can be tuned from topologically nontrivial states to topologically trivial states by varying the orientation angle Φ.

Under C4 rotational symmetry, the lattices remain invariant upon a 90° rotation. Accordingly, three orientation angles Φ=0°, 45°, and 70° are selected for analysis. Figures 3(a)–3(c) show the distributions of Re(Ez), where the black arrows represent the directions of the in-plane electric-field vectors. The symmetry of lattices is further confirmed by the Fourier transform of Re(Ez), shown in the insets of Figs. 3(a)–3(c). The electric-field meron lattices in Figs. 3(a) and 3(b) exhibit a Néel-type topology at Φ=0° and a twisted-type topology at Φ=45°, respectively. However, at Φ=70°, the calculated skyrmion number is 0.365, and the electric-field vectors exhibit a discontinuous distribution, indicating a nontopological lattice, as shown in Fig. 3(c). The mechanism underlying the phase transitions in the electric-field meron lattice can be understood through the propagation constant diagram in Fig. 1(b). Similar to the case of TM SPP modes at the metal/isotropic media interface, the surface waves at the metal/uniaxial crystal interface remain a pure TM SPP mode at Φ=0°, 90°, and 180°, which supports the formation of the electric-field meron lattices. However, at other orientation angles, the surface waves evolve into hybrid SPP modes, leading to the observed topological phase transitions. Notably, at Φ=45° and 135°, the off-diagonal components of the dielectric tensor of the uniaxial crystal reach their maximum values, resulting in strong coupling between TE and TM SPP modes. This strong coupling leads to the propagation constant β reaching a peak and gives rise to the twisted-type topological lattice. Figures 3(d)–3(f) represent the corresponding topological textures of electric-field vectors encoded in the HSL color space: saturation is maximal, lightness encodes the nz component, and hue corresponds to arctan(ny/nx). These visualizations reveal that the electric-field lattice becomes structurally unstable as Φ varies, revealing that the OA orientation governs the topological phase transitions of the electric-field meron lattices.

Then, we consider the formation of topological lattices for the quantum number of TAM l=2. It is worth noting that, at the metal/isotropic media interface, the topological lattice formed at l=2 is similar to that of l=0, as the propagation constants of all OA orientation angles are the same. However, at the metal/anisotropic media interface, the situation is different. At Φ=0°, an electric-field meron lattice is formed, with its 3D vector distribution depicted in Fig. 4(a). It can be found that the topological lattice formed at l=2 differs from that of l=0. This difference arises from the distinct quantum number in phase terms of Eq. (4): Ez∝cos(β1y)+cos(β2x) for l=0 and Ez∝cos(β2x)−cos(β1y) for l=2. Furthermore, we calculate the skyrmion number of the marked unit cell via Φ [Fig. 4(b)]. A significant transition occurs near Φ=45°, where the skyrmion number undergoes a sharp gradient. This sharp variation demonstrates strong sensitivity to the OA orientation, which is potentially applied in the optical sensing with ultrahigh sensitivity. Figures 4(c) and 4(d) illustrate the distributions of Re(Ez) and their corresponding topological textures for Φ=0°, 45°, and 90°, respectively. The calculated skyrmion numbers are Q=−0.490, 0, and 0.490, respectively. Notably, comparing Fig. 4(c3) with Fig. 4(c1), a reversal of the electric-field vector directions is observed within the marked unit cell at Φ=90°. This reversal originates from the difference in electric-field vector orientation induced by Φ in phase term exp(ilϑn) of Eq. (4), and explains why the skyrmion number becomes positive in the range of Φ approximately 45°–135°, as illustrated in Fig. 4(b). It can be observed that it even caused the annihilation of topological lattice due to the appearance of an imaginary unit in the electric-field component Ez∝i[cos(2β2x/2)−cos(2β1y/2)] when Φ=45° at t=2π/ω. It should be emphasized that for l=2, when the OA orientation angle is not an integer multiple of 0° or 45°, the skyrmion number calculated from n=Re(E)/|E| is not strictly accurate. In such cases, the electric-field components are not in phase, and due to the presence of temporal components, the field oscillates at different time points. As a result, the electric-field structure varies with time, making it impossible to form a strictly defined skyrmion. However, since only the real part of an electric-field component can be measured experimentally, we use n=Re(E)/|Re(E)| to calculate the skyrmion number here.

Finally, we consider the case of nonzero quantum number of TAM. In evanescent vortex fields, the formation of spin lattices induced by SOC has already been reported [20,21]. As for the surface plane wave at the metal/uniaxial crystal, the time-average Poynting vector is given by P=Re(E*×H)/2, where the asterisk denotes the complex conjugate, whereas the kinetic momentum of light associated with the time-average Poynting vector is p=Re(D*×B)/2. Then, the TAM is J=r×p with r the position vector, and it can be re-expressed by the sum of SAM and orbital angular momentum (OAM) through vector calculation with identities based on Eq. (B2) in Appendix B as [64] $$\begin{gathered} ∫−∞∞Jd3x=∫−∞∞r×pd3x \\ =∫−∞∞{r×14ωIm[B*·(∇)H+D*·(∇)E]+14ωIm(B*×H+D*×E)}d3x. \end{gathered}$$(6)

The first term in the right-hand side of Eq. (6) is the OAM density, and the second term is the SAM density.

We perform a comparative analysis of the cases with the quantum numbers of TAM l=1 and l=3. Under C4 rotational symmetric field, the SOC present in the hybrid surface waves field leads to the formation of spin meron lattices by the unit vectors of SAM as n=S/|S|.

In both cases, Figs. 5(a) and 5(b) depict the normalized out-of-plane SAM distributions Sz at Φ=0°, which exhibit C4 symmetry. The corresponding 3D spin orientation distributions show the formation of spin meron lattices, as shown in Figs. 5(c) and 5(d). Within each unit cell, the local spin orientations gradually transition from a central “up” or “down” state to an edge flat state (Sz=0). The calculated skyrmion numbers within the black-framed regions are Q=0.498 and Q=−0.498, respectively. To evaluate the stability of these lattices, we calculate the skyrmion number of the central unit cell under different Φ, as shown in Fig. 5(e). It can be observed that regardless of Φ changes, the skyrmion numbers remain nearly constant at approximately ±0.5 in both cases. This remarkable stability indicates that the topological spin lattice exhibits robustness to varying Φ.

To further explore the topological characteristics of spin lattices for l=1 and l=3, we chose Φ=0°, 45°, and 70°, as the instance of l=0. The first row of Figs. 6(a) and 6(b) presents Sz distributions, with black arrows indicating the directions of the normalized in-plane SAM. The insets show the Fourier transforms of Sz to observe the symmetry of lattices. The second row illustrates the calculated phase distributions of Ez, where black arrows represent the energy flow of the hybrid surface waves, exhibiting alternating circulations. Specifically, each vortex center of the Poynting vector corresponds to a phase singularity of Ez. The third row visualizes the spin vector orientations, clearly displaying the evolution of topological spin textures. At Φ=0° and 45°, each unit cell exhibits a Néel-type meron texture, whereas at Φ=70°, a twisted-type meron texture emerges. Despite these changes in spin texture, the skyrmion number remains unchanged, reflecting its topological protection. This demonstrates both the invariance of the skyrmion number and the robustness of spin meron lattice textures against changes in Φ.

The origin of this topological stability can be attributed to two key factors. First, as shown in Eq. (6), the decomposition of TAM into orbital and spin components in an anisotropic medium is similar to that in an isotropic medium. Second, the results illustrated in Fig. 6 indicate that within a uniaxial crystal, the directions of spin and energy flow density are locked with each other to maintain the stability of the skyrmion number of topological textures. According to the spin-momentum equation S∝∇×p [62], the surface integral of SAM over the transverse plane can be converted into a line integral of the kinetic momentum p. As the integration region extends to infinity, p approaches zero, so the surface integral of S is also zero. Consequently, SAM must simultaneously exhibit positive and negative values in a continuous distribution. This ensures the existence of regions where spin vector orientation transitions from positive to zero and from negative to zero, thus maintaining the topological stability of the spin meron lattice.

In conclusion, we have demonstrated the formation of optical meron lattices at the metal/uniaxial crystal interface. Specifically, electric-field meron lattices are formed when the TAM quantum numbers are l=0 and l=2, while spin meron lattices emerge for l=1 and l=3, and the topological lattices exhibit periodicity (=4) with respect to the quantum number of TAM. Notably, the skyrmion number of the spin meron lattices is invariant for various values of Φ, indicating its topological stability. This invariance can be attributed to the fact that the decomposition of TAM into orbital and spin components in an anisotropic medium is similar to that in an isotropic medium. Furthermore, based on the spin-momentum relations, the SAM exhibits a continuous distribution containing both positive and negative values, thereby ensuring the topological stability of the spin meron lattice. In contrast, as the orientation angle Φ is tuned, the surface waves are no longer pure TM SPP modes at the metal/uniaxial crystal interface, leading to topological phase transitions in the electric-field meron lattices. These results provide novel physical insights into the topological stability and phase transitions of topological quasiparticles in optical systems, and offer an avenue for exploring applications in advanced spin-optics and topological photonics.