Designing rotational motion of charge densities on plasmonic nanostructures excited by circularly polarized light

2.1 C ₄ symmetry

Model calculations were performed to predict the rotational motion of charge densities in plasmonic nanostructures with different rotational symmetries. First, we consider a plasmonic system with C ₄ symmetry and its azimuthal periodic plasmon modes. The charge density for the azimuthal periodic plasmon mode excited by x-polarized linearly polarized (LP) light, ρ m x p o l ( r , θ ) , is expressed in polar coordinates using radius r and angle θ as follows.

where A(r) denotes the radial distribution. The variable m is a positive integer that represents the order of the rotational symmetry of the charge density, where m = 1, 2, and 3 correspond to the dipole, quadrupole, and hexapole modes, respectively [20].

In a system with C ₄ symmetry, the charge density excited by y-polarized LP light, ρ m y p o l ( r , θ ) , is required to be exactly equal to ρ m x p o l rotated by π/2. This can be expressed as

By classifying the odd-order resonance modes into m = (4j ± 1) using integer j, we can divide Eq. (2) as follows:

Note that even-order modes are forbidden in C ₄ symmetric systems. The modeled ρ m x p o l and ρ m y p o l for the dipole mode, m = 1 = 4 ⋅ 0 + 1 (j = 0), and those of the hexapole mode, m = 3 = 4 ⋅ 1 − 1 (j = 1), are plotted in Figure 2(a) and (b), respectively. A(r) is set to 1 for simplicity.

Figure 2:

Model calculations of charge density distributions for C ₄ symmetry. (a and b) Simplified charge density distribution of the (a) dipole and (b) hexapole mode excited by x-and y-polarized linearly polarized (LP) light, respectively. The black arrows indicate the direction of the electric field of excitation light. (c and d) Charge density distributions with ϕ = π/4, calculated from the charge density of each of the x- and y-polarized excitation. The yellow circular markers and purple arrows indicate the angular shift where the charge density is maximum. (e and f) Angular distribution of charge density under circularly polarized excitation for one oscillation period.

The charge density distribution excited by the LP light with a polarization angle ϕ, ρ _m (r, θ; ϕ), is expressed as a superposition of the charge density distribution excited by x-polarized light and that excited by y-polarized light,

From Eqs. (1)–(4), the charge density distribution of order m = 4j ± 1 can be transformed as follows.

Equation (5) indicates that the charge density distribution ρ _m (r, θ; ϕ) is equal to ρ m x p o l ( r , θ ) rotated by ϕ/m counterclockwise (clockwise) for m = 4j + 1 (m = 4j − 1).

In Figure 2(c) and (d), we plot the calculated charge density distribution excited by LP light with polarization angle ϕ = π/4, ρ ₁(θ; ϕ = π/4) and ρ ₃(θ; ϕ = π/4). In the dipole mode plotted in Figure 2(c), the angle at which the charge density exhibits its maximum value is oriented π/4, aligning with the direction of the electric field. Conversely, in the hexapole mode, the maximum charge density shifts in the opposite direction to the shift of the electric field by −π/12.

It is noteworthy that the charge density beneath the black arrow (aligned with the electric field vector) is positive for both horizontal and vertical polarizations (Figure 2(b)), but turns negative for diagonal polarization (Figure 2(d)). This non-intuitive charge density distribution cannot occur in systems with C ₈ symmetry, such as regular octagons and isotropic disks, where ρ(θ − π/4) must be ρ(θ) rotated clockwise by π/4. C ₄-symmetric systems lack C ₈ symmetry, resulting in horizontal and diagonal excitations producing fundamentally distinct outcomes (i.e., ρ(θ − π/4) differs from ρ(θ) rotated clockwise by π/4).

Under CP light irradiation, ϕ continuously varies from 0 to ±2π within an oscillation period. Figure 2(e) and (f) show the time evolution of the angular charge density distribution for modes m = 1 and 3, over one oscillation period (ϕ varies continuously from 0 to 2π), corresponding to the incidence of left-handed CP light. For the dipole mode, the angle at which the charge density attains its maximum value, indicated by the yellow marker and purple arrow, follows the rotation of the electric field and rotates 2π per oscillation period. In contrast, the angle of maximum charge density in the hexapole mode rotates in the direction opposite to the rotation of the electric field, with a rotation angle of 2π/3 per oscillation period.

2.2 C _N symmetric system

We now extend our discussion to general rotationally symmetric systems. The C _N rotationally symmetric systems require the following symmetry, using ϕ _N = 2π/N,

where, ρ m N ( r , θ ; ϕ N ) represents the charge density excited by the LP light with polarization angle ϕ _N . From Eq. (4), ρ m y p o l ( r , θ ) can be expressed as

By substituting m = (Nj ± 1), the first term of the numerator in Eq. (7), cos[m(θ − ϕ _N )], can be transformed to

because Njϕ _N = 2πj. From Eqs. (7) and (8), following expression can be derived for the case m = Nj ± 1, as in Eq. (3),

Here, ρ m x p o l , ρ m N , and calculated ρ m y p o l for the case of m = N ± 1 (j = 1) in the system with rotational symmetry of N = 3 (m = 2, 4) and N = 6 (m = 5, 7) are plotted in Figure 3(a). Note that the triangular and hexagonal frames here, and the square frames in Figure 2, are introduced as typical representations of each rotational symmetry and do not indicate a specific geometry. Owing to Eqs. (3) and (9) being identical except for the N, the angular distribution of the m-th order mode excited by the LP light with angle ϕ for the C _N rotational symmetry can also be expressed as

Figure 3:

Model calculations of charge density distributions for C _N symmetry. (a) Calculated charge density distributions in the C ₃ (top panel) and C ₆ (bottom panel) symmetric system. From the left to right column, the polarization of the excitation lights is ϕ = 0, ϕ = ϕ _N = 2π/N, and ϕ = π/2. (b) Calculated charge density distributions for ϕ = π/4. The yellow markers and purple arrows indicate the angular shift where the charge density is maximum.

Considering the rotation of the electric field for CP light, Eq. (10) indicates that the rotational charge motion in a generalized rotational system can be described as follows: the charge density distribution of azimuthal periodic modes with m = Nj ± 1-order excited in C _N symmetric system rotates by ±2π/m in one oscillation period.

Calculated charge density distributions, ρ _m (θ; ϕ = π/4) for the m = N ± 1-th modes, are plotted in Figure 3(b). Similar to the C ₄ system shown in Figure 2, the rotational direction of the charge density distribution for each mode is aligned with that of the electric field for m = Nj + 1 and opposite for m = Nj − 1. The amount of the rotation angle is ϕ/m in both cases.

It should be noted that the direction and speed of rotation of the charge density obtained here are consistent with one of the accessible OAM values of the far-field scattered field from a rotationally symmetric scatterer excited by CP light [26]. The relationship between the rotational symmetry of the scatterer and the possible OAM values of the far-field scattered field can be understood based on the concept of angular grating, which extends the linear grating to azimuthal periodic polar systems [15], [27], [28], [29]. Our calculations confirm the similarity between the rotational motion of charge densities within the plasmonic structure and far-field scattered field, both governed by the same restriction of rotational symmetry.

3.1 Isolated nanoantenna

To confirm the validity of the model calculations, we evaluated the temporal evolution of the near-field electromagnetic distribution around plasmonic structures using FDTD simulations. Figure 4(a) shows a schematic of the model used in the FDTD simulation. A 50 nm thick Au plate with side length L = 1 μm is placed in the xy-plane, and CP light with wavelength λ _inc is normally incident from a total-field scattered-field (TFSF) source located 5 μm away. The boundary conditions for the x, y, and z-directions are perfect matched layers (PML), and the minimum mesh width is 10 nm. The dielectric function of Au is obtained from the Palik database [30]. The observation plane is set in the middle of the Au layer to record electric field E and the magnetic field H . All simulations were performed using ANSYS Lumerical (2023 R1). The charge density distribution generated by a surface plasmon wave propagating along a planar surface is well-characterized by the plasmon’s electric and magnetic field amplitudes in the direction perpendicular to the propagation [31], [32]. The electric and magnetic amplitudes reach their peak values where the charge density is highest. Keeping this in mind, we examine the H _z distribution for our plasmon waves orbiting in the xy plane (as shown in Figure 4(a)). We see that the spatial distribution of H _z component corresponds well with that of the charge density orbiting within the plasmonic structures (see Supplementary Material).

Figure 4:

FDTD simulations for a square nanoantenna. (a) Schematic of the FDTD simulation model showing the setup with a wave source, an Au plate in the xy-plane, and an observation plane. (b) Spatial distributions of the H _z component at a representative phase for a square plate with L = 1 μm, excited by circularly polarized (CP) light. (c) Spatial phase distributions of the H _z component plotted in (b).

The spatial distributions of the H _z components at the characteristic wavelengths in the square nanoantenna are plotted in Figure 4(b). From these distributions, we determined that the resonance orders m are 1, 3, 5, and 7 at wavelengths of 3 μm, 1.3 μm, 0.89 μm, and 0.68 μm, respectively. Trajectories of the positive H _z components during one oscillation period, evaluated from the time evolution of the magnetic field distribution for each wavelength, are plotted by black arrows in Figure 4(b). (Supplementary Movie 1 shows an animation of |E| and H _z components for one oscillation period). Their direction and length indicate that the rotational directions at adjacent resonance orders are opposite each other, and that the amount of rotation per oscillation period is 2π/m. These direction and speed of the rotation for each resonance mode are consistent with those obtained from Eq. (10).

Analogous to how the sign and magnitude of the OAMs of a vortex beam can be deduced from the spatial phase distribution of the electromagnetic field, the direction and speed of the rotational motion of the charge density on plasmonic structures can be evaluated from the spatial phase. Figure 4(c) shows the spatial phase distribution of the H _z component for each wavelength plotted in Figure 4(b). By examining the phase gradient along the perimeter of the square structure, both the direction of the phase gradient and the number of phase cycles per perimeter are observed to be dependent on the resonance order. These continuous, unidirectional phase gradients from −π to π indicate that the charge density that generate the H _z components are in continuous orbital motion around the perimeter of the structure.

The phase gradients along the perimeter are represented by red arrows for clockwise and blue arrows for counterclockwise. The direction of the phase gradient indicate that the rotational direction of the charge densities in the m = 4j + 1 (m = 4j − 1) mode is the same (opposite) as that of the handedness of the excitation CP light. The region indicated by the dashed line, where the regularity of the phase gradient is disturbed, is attributable to the interference from adjacent off-resonant modes with peaks at close frequencies. In particular, at high frequencies, where the resonance frequencies are close to each other, the electromagnetic field distribution at each wavelength is affected by adjacent modes rotating in opposite directions.

3.2 Rectangular nanoantenna

In C _N systems where N ≥ 3, m is uniquely represented by N and j, and thus the rotational direction is automatically determined from Eq. (10) for each resonance order. However, for a C ₂ symmetric system, there are two ways to represent m, and the direction of the rotational motion of the charge density is not determined solely by rotational symmetry. For instance, the m = 1 case in the C ₂ symmetric system can be represented by m = 1 = 2 ⋅ 0 + 1 for j = 0 and m = 1 = 2 ⋅ 1 − 1 for j = 1. Figure 5(a) shows modeled ρ 1 x p o l and ρ 1 N ( ϕ = π ) in the system with rotational symmetry of N = 2. The pairs of ρ defined by rotational symmetry do not form angles and consequently cannot serve as a basis for ρ(ϕ) excited at arbitrary polarization angles ϕ. Instead, ρ ^ypol is determined independently of ρ ^xpol in C ₂ symmetric systems. Using Eq. (3), the charge density ρ m y p o l ( r , θ ) for the y-polarized excitation in the above instance is estimated to be +A(r) sinθ for j = 0 and −A(r) sinθ for j = 1. Consequently, the charge density can rotate either to the clockwise or counter-clockwise, as shown in Figure 5(b)

Figure 5:

Model calculations and FDTD simulations for a rectangular nanoantenna. (a and b) Modeled charge density distributions in the C ₂ symmetric system. From the left to right column, the polarization angle of the excitation light is (a) ϕ = 0, π, (b) π/2, and π/4. The yellow markers and purple arrows indicate the angular shift where the charge density is maximum. (c) Spatial distributions over half an oscillation period and spatial phase distributions of H _z components for a rectangular nanoantenna at λ = 5 μm (top panels) and 2.5 μm (bottom panels).

In practice, such conditions can be designed by the geometry of the structure, or by the phase delay of the plasmonic charge motion with respect to the external field. In the case of rectangular nanoantennas, the long and short axes, having different resonance wavelengths, λ _lo and λ _sh, exhibit distinct phase delays relative to the external field [33]. When the external field wavelength (λ _inc) is much longer than both λ _lo and λ _sh, the plasmonic resonances on both axes are in phase with the external field. Conversely, if λ _inc lies between λ _lo and λ _sh, the resonances on the short axis oscillate in phase, whereas those on the long axis oscillate out of phase.

The H _z components of a rectangular nanoantenna excited by CP light are plotted in Figure 4(d), over a half oscillation period. The initial phase (ϕ ₁, ϕ ₂) is defined where the positive H _z is localized on the left side of the rectangle at each wavelength. The trajectories of positive H _z components for the half oscillation period are plotted with black arrows in the fourth image. The lengths of the long and short axes of the rectangle are 1.2 and 0.4 μm, respectively, and the evaluated resonance wavelengths are λ _lo = 3.6 μm and λ _sh = 1.8 μm. At λ _inc = 5 μm, which is a longer wavelength than both λ _lo and λ _sh, the H _z distribution rotates clockwise along the perimeter of the rectangle, which is the same as the handedness of the CP light (top panels). In contrast, at λ _inc = 2.5 μm, which is located between λ _lo and λ _sh, the direction is reversed and counterclockwise rotation is observed (bottom panels). The phase distribution of H _z shown in right panels of Figure 4(d) clearly indicates phase gradients in opposite directions.

These results confirm that the rotational behavior of charge densities in isolated plasmonic nanoantennas is determined by the rotational symmetries of the structure and the resonance mode, as calculated using Eq. (10). In particular, the results for the rectangular nanoantenna exhibit that the geometric design of the plasmonic structure allows for the selective excitation of the rotational motion of the charge density.

3.3 Periodic lattice

As the only requirement for Eq. (10) is rotational symmetry, the discussion of the rotational behavior of charge densities is applicable to systems other than isolated plates, as long as they are rotationally symmetric. Circular structures, which generally allow only dipole mode for normally incident plane wave excitation because of their infinite rotational symmetry [19], [34], can exhibit multipole modes when their symmetry is restricted by an array structure [35], [36]. In other words, the multipole modes excited in the periodically arranged structure arise from the external rotational symmetry introduced by the arrangement. Therefore, the multipole modes excited by CP light in these periodic structures should also exhibit rotational behavior determined by the rotational symmetries.

Here we introduced C ₄ and C ₃ symmetry by arranging circular apertures in square and hexagonal lattice arrays, respectively, as illustrated in Figure 6(a). As with the isolated plate, the resonance order and rotational direction of the excited plasmon mode can be evaluated from the H _z component on the inner wall surface. From the simulation model for the isolated plate used in the previous section, the boundary condition in the xy direction is changed to a periodic boundary, and the TFSF wave source is replaced with a plane wave source. The lattice structure comprises a 50 nm thick Au thin film in vacuum, featuring an array of holes with radius r and lattice constant d. The minimum mesh size is reduced to 5 nm to minimize the effect of unintended C ₄ symmetry introduced by the Yee lattice in the FDTD simulation.

Figure 6:

FDTD simulations for hole array structures. (a) Schematic of the lattice structures used in the FDTD simulations. (b–d) Spatial distribution of the (b) H _z components at representative phase, (c) phase distributions of the H _z components, and (d) Poynting vectors P (blue arrows), time-averaged current densities J _dc (red arrows), and averaged intensity distributions of the electric field | E | (color map). Calculations were performed with r = 350 nm and d = 800 nm for square lattices, and r = 400 nm and d = 1,200 nm for hexagonal lattices. The scale bars in (b) represent 200 nm.

Figure 6(b) shows the spatial distributions of the H _z components at representative phases for different resonance orders. The trajectories of the positive H _z components per oscillation cycle are also plotted by black arrows, as in Figure 4(b). (Animation of | E | and H _z components for one oscillation period are shown in Supplementary Movies 2 and 3). The spatial phase distributions of H _z components are plotted in Figure 6(c), and the phase gradient along the inner wall of the circular aperture is indicated by red and blue arrows. The rotation angle indicated by each arrow is approximately 2π/m and is opposite to the incident CP light only when m = Nj − 1. This observation is consistent with Eq. (10).

It is noteworthy that the m = 5 mode corresponds to N + 1 for N = 4 and 2N − 1 for N = 3, indicating that multipole modes of the same order rotate in opposite directions on square and hexagonal lattices. These calculations demonstrate that the rotational symmetry introduced externally by the periodic array allows us to manipulate the uniform rotational motion of charge densities in a wide range of plasmonic structures with a high degree of freedom using spatially homogeneous CP light.

One might think that the observed rotational motion of the charge densities could be due to the wagon wheel effect resulting from the discrete computation time. Additionally, the appearance of a rotating electromagnetic field does not explicitly indicate that the plasmon wave is in rotational motion with unidirectional energy flow, as seen in instances where the group velocity of a phase-adjusted structured light can propagate opposite to its energy flow direction [37]. To demonstrate more explicitly the direction of energy flow, we calculate the Poynting vector and current densities.

The time-averaged Poynting vector P , representing the energy flow in an electromagnetic field, and the time-averaged current density J _dc, representing the charge flow within the metal, were calculated using the following equations, respectively [38], [39]. The derivation of J _dc is summarized in the Supplementary Material.

where e is the charge of the electron, n is the charge density at rest, and σ _ω is the dynamic conductivity of the metal.

Figure 6(d) shows the in-plane components of the Poynting vectors (blue) and current densities (red), plotted on the color map of the averaged electric field intensity. The anisotropic intensity distributions are attributed to interference from adjacent off-resonant modes (see Supplementary Material). The orientation of the Poynting vectors and time-averaged current densities along the circumference corroborates that the observed dynamics in the electromagnetic field distribution is reflecting the continuous orbital motion of charge density or a surface current.

Manipulating the rotational motion of charge density distributions in periodic plasmonic nanostructures shows that uniform control of rotational charge motion in large-area structures suitable for experimental observation is achievable. This capability offers a variety of potential applications, such as the observation of the inverse Faraday effect [16], [17], [18], [39], [40], designing the resonant interaction between magnetic resonant antennas [41], manipulation of the photocurrent [42], [43], [44] and spin currents [45], [46], [47], [48], [49], [50], [51], and control of the SAM associated with the plasmonic near-field [31], [52], [53], [54].

The spatial phase distributions similar to those of vortex beams with higher-order OAMs observed in Figures 4(c) and 6(c) suggest the possibility of experimental investigations of light–matter interactions involving OAM [55], [56], [57], [58], [59]. Furthermore, our calculations indicate that the rotational charge motions on plasmonic structures excited by CP light exhibit opposite directions between adjacent modes among the multipole resonance modes that can be excited in each structure. This is important for explaining certain complex spatial distributions and non-intuitive temporal variations of electromagnetic fields often observed in electromagnetic simulations using CP excitation light [60], and aids in the calculation of the optical properties of enhanced electromagnetic fields, such as SAM [52] and optical chirality [61].

Although we limited the structures in each calculation to one type for simplicity, more complex rotational behavior of the charge densities can be designed by combining different structures. For instance, charge densities can be rotated in opposite directions on the same plane by using structures of different sizes that share the same rotational symmetry (see Supplementary Movie 4). Designing geometries other than regular polygons, such as swastikas, would facilitate a more efficient excitation of specific rotational behaviors. The actual charge density distribution in nanopolygon plates often does not exhibit the perfect azimuthal periodicity as simplified in this study [62], [63]. Defining basis vectors based on modal analysis that considers the geometry of each structure to evaluate such complex distributions is our future work [64].