Optical computation of divergence operation for vector field

Spatially varying polarization is one of the most fundamental properties of optical vector fields [1]. Recently, more and more studies show that the polarization singularities in vector fields play the key role to understand topological physical phenomena. Especially, in a periodic dielectric slab, a vortex as the winding of far-field polarization vectors in k space connects to the bound state in the continuum (BIC) and forms ultra-high-Q resonances in open systems [2, 3, 4, 5, 6, 7], which have been exploited for BIC laser [8]. Also for electromagnetic multipoles it is important to resolve such singularities in the polarization states for identifying the closed radiation channels [9].

In order to analyze such vector fields, spatial polarization distributions are commonly measured by mechanically rotating polarizers and wave plates, which is time-consuming and unstable to infer the polarization singularities from the intensity variation. Over the past few years optical analog computing has attracted particular attentions with the advantage of real-time and high-throughput computations [10]. Various optical analog computing devices were designed for spatial differentiation for scalar fields simply with linear polarization, e.g. metamaterials/metasurface [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], dielectric interface or periodic slabs [22, 23, 24, 25, 26, 27, 28, 29, 30, 31], surface plasmonic structures [32, 33, 34, 35]. Therefore, it can be intriguing to extend the optical computing method from scalar operation to vector one, which is capable of both directly analyzing the divergence of vector fields and revealing the polarization singularities, on a single shot.

In this article, a periodic plasmonic metasurface is proposed to perform optical analog computation of divergence operation. For linearly polarized light, our previous studies have experimentally demonstrated the spatial differentiation by two different approaches. First, based on the simplest planar plasmonic structure, the spatial differentiation along the direction parallel to the incident plane can be realized by spatial mode interference when surface plasmon polariton (SPP) is excited [32, 34, 33]. Second, with the spin Hall effect of light, by analyzing polarization states orthogonal to incident light, the spatial differentiation perpendicular to the incident plane can be realized during reflection at a single optical planar interface, which generally occurs at any planar interface, regardless of material compositions or incident angles [29, 30]. Here in order to realize optical computing of divergence operation for vectorial paraxial beams, we introduce periodic slots in the metasurface and combine these two effects together.

Such a metasurface simultaneously computes spatial differentiation along the directions parallel and perpendicular to the incident plane. We enable the divergence computation by creating the constructive interference between two derivative results. Generally, for vectorial paraxial beams, the electric fields are dominant in transverse directions. Since the divergence of total electric fields are always absent for the light propagating in air, the divergence of vectorial paraxial fields is contributed by the nonzero derivative of the longitudinal components along the beam propagation direction, which corresponds to the polarization singularities of sources or sinks.

The periodic plasmonic metasurface is depicted in Fig. 1(a), where a silicon slot grating is fabricated on a thick gold layer substrate and then a monolayer graphene is transferred onto the grating. Here the air-graphene-silicon interface supports SPP modes which are strongly confined and propagate along the graphene layer in THz frequency region [34]. The SPP on graphene is excited through phase matching enabled by silicon grating: ${k_{0}}\sin\theta-{{2\pi}\mathord{\left/{\vphantom{{2\pi}\Lambda}}\right.\kern-1.2pt}\Lambda}=-{\beta_{\rm{spp}}}$, where $k_{0}$ is the wavevector in vacuum, $\theta$ is the incident angle, ${\beta_{\rm{spp}}}$ is the wavevector of SPP. We note that here the phase matching condition ensures that only the zeroth order beam is diffracted.

We now consider that a vectorial paraxial beam illuminates on the metasurface with the incident electric fields ${{\bf{E}}_{\bf{i}}}=E_{x}^{i}\left({x,y}\right){{\bf{e}}_{\bf{x}}}+E_{y}^{i}\left({x,y}\right){{\bf{e}}_{\bf{y}}}$, where $x$ and $y$ are the beam coordinates perpendicular to the propagation direction, as depicted in Fig. 1. Correspondingly, the reflected beams are written as ${{\bf{E}}_{\bf{r}}}=E_{x}^{r}\left({x,y}\right){{\bf{e}}_{\bf{x}}}+E_{y}^{r}\left({x,y}\right){{\bf{e}}_{\bf{y}}}$. Here $E_{x}^{i(r)}$ and $E_{y}^{i(r)}$ represent $x$ and $y$ components of the vectorial paraxial fields, for the incident (reflected) beams, respectively. With the convention of $\exp\left({-i\omega t}\right)$, each component can be decomposed to a series of plane waves: $E_{\alpha}^{i\left(r\right)}\left({x,y}\right)=\iint{\tilde{E}_{\alpha}^{i\left(r\right)}\left({{k_{x}},{k_{y}}}\right){e^{i{k_{x}}x}}{e^{i{k_{y}}y}}d{k_{x}}d{k_{y}}}$. Since the diffracted light of the metasurface is only in the zeroth order, each plane wave with the wavevector component $\left({{k_{x}},{k_{y}}}\right)$ only generates the reflected plane wave with the same $\left({{k_{x}},{k_{y}}}\right)$. Therefore, the reflected vector fields are determined by four spatial spectral transfer functions ${H_{\alpha\beta}}\left({{k_{x}},{k_{y}}}\right)$ defined as ${H_{\alpha\beta}}\left({{k_{x}},{k_{y}}}\right)={{\tilde{E}_{\alpha}^{r}\left({{k_{x}},{k_{y}}}\right)}\mathord{\left/{\vphantom{{\tilde{E}_{\alpha}^{r}\left({{k_{x}},{k_{y}}}\right)}{\tilde{E}_{\beta}^{i}\left({{k_{x}},{k_{y}}}\right)}}}\right.\kern-1.2pt}{\tilde{E}_{\beta}^{i}\left({{k_{x}},{k_{y}}}\right)}}$, where $\alpha$, $\beta$ represent $x$ or $y$ electrical components.

Next we show that the spatial spectral transfer functions ${H_{\alpha\beta}}\left({{k_{x}},{k_{y}}}\right)$ can be engineered to realize the optical computing of divergence operation. First, when the wavevector of an incident plane wave is out of the incident plane, i.e. ${k_{y}}\neq 0$, the incident $y$-component field can excite the $x$-component reflected field, and vice versa, due to the transversal mode requirement of plane waves [29, 30]. Especially when the grating slot is shallow, that is, the grating is considered as a perturbation to the air-graphene-silicon interface, under the paraxial approximation, the spatial spectral transfer function ${H_{xy}}$ has the form of

$${H_{xy}}=i{A_{y}}{k_{y}}$$

(1)

where ${A_{y}}$ is a complex constant. Eq. (1) indicates that when the incident paraxial beam is dominant along the $y$ direction ${{\bf{E}}_{\bf{i}}}=E_{y}^{i}\left({x,y}\right){{\bf{e}}_{\bf{y}}}$, by analyzing the reflected field along the $x$ direction, the output polarized beam is $E_{x}^{r}={A_{y}}\frac{{\partial E_{y}^{i}\left({x,y}\right)}}{{\partial y}}$, and ${A_{y}}$ represents a scaling factor for the $y$-directional differentiation operation.

On the other hand, through phase matching, the incident electric field component in the incident plane can excite the SPP wave at the air-graphene-silicon interface. Meanwhile the excited SPP field propagating along the surface of grating leaks out and generates radiation field through grating coupling process. The total reflected field in the incident plane is contributed by the interference of direct reflected wave and leakage of SPP mode. Such a coupling process can be described by the spatial coupled-mode theory [32, 36]. Especially, when the loss rate due to the SPP leaky radiation is equal to that of the material absorption, i.e. under the critical coupling condition, the total reflected field is the first-order spatial differentiation of the incident one along the $x$ direction [34]. In this situation, the corresponding spatial spectral transfer function is

$${H_{xx}}=i{A_{x}}{k_{x}}$$

(2)

where ${A_{x}}$ is a complex constant determined by the loss rate due to the SPP leaky radiation. Eq. (2) shows that when the incident paraxial beam is dominant along the $x$ direction ${{\bf{E}}_{\bf{i}}}=E_{x}^{i}\left({x,y}\right){{\bf{e}}_{\bf{x}}}$, by analyzing the reflected field along the $x$ direction, the output polarized beam is $E_{x}^{r}={A_{x}}\frac{{\partial E_{x}^{i}\left({x,y}\right)}}{{\partial x}}$, where ${A_{x}}$ represents the scaling factor for the $x$-directional differentiation operation.

Therefore, for an incident vector field ${{\bf{E}}_{\bf{i}}}=E_{x}^{i}\left({x,y}\right){{\bf{e}}_{\bf{x}}}+E_{y}^{i}\left({x,y}\right){{\bf{e}}_{\bf{y}}}$, after analyzing the reflected field along the $x$ direction, the output polarized beam $E_{x}^{r}$ is the interference result between the derivatives of two incident field components $E_{x}^{i}\left({x,y}\right)$ and $E_{y}^{i}\left({x,y}\right)$ along the $x$ and $y$ directions respectively. We note that generally the scaling factors ${A_{y}}$ and ${A_{x}}$ are two complex numbers with different phases. Suppose the wave plate in Fig. 1(a) has a phase delay $\Delta\phi$ in the $y$ component, the output polarized beam after analyzing the reflected fields along the $x$ direction is

$$E_{x}^{r}\left({x,y}\right)={A_{x}}\frac{{\partial E_{x}^{i}\left({x,y}\right)}}{{\partial x}}+{A_{y}}{e^{i\Delta\phi}}\frac{{\partial E_{y}^{i}\left({x,y}\right)}}{{\partial y}}$$

(3)

Furthermore, when a constructive interference between two derivatives with the same scaling factor is realized, i.e. $\left|{{A_{x}}}\right|=\left|{{A_{y}}}\right|$ and $\Delta\phi=\arg\left({{A_{x}}}\right)-\arg\left({{A_{y}}}\right)$, we enable a vector-to-scalar field transformation which corresponds to the divergence computation of the incident vector fields as

$$E_{x}^{r}\left({x,y}\right)={A_{x}}{\nabla_{\bot}}\cdot{{\bf{E}}_{{\bf{i}}}}$$

(4)

where ${\nabla_{\bot}}={{\bf{e}}_{\bf{x}}}\frac{\partial}{{\partial x}}+{{\bf{e}}_{\bf{y}}}\frac{\partial}{{\partial y}}$ is the transverse divergence operator.

To demonstrate the divergence computation, we numerically simulate such a vector-to-scalar field transformation by the finite element method using the commercial software COMSOL. Here the surface conductivity of monolayer graphene is described by the Drude model $\sigma=\frac{{{e^{2}}}}{{\pi{\hbar^{2}}}}\frac{{i{E_{\rm{F}}}}}{{\omega+i{\tau^{-1}}}}$, where $e$ is the elementary charge, $\tau$ is the carrier relaxation time which satisfies $\tau=\mu{E_{\rm{F}}}/ev_{\rm{F}}^{2}$. ${E_{\rm{F}}}$, $\mu$ and ${v_{\rm{F}}}$ correspond to the Fermi energy, the carrier mobility and the Fermi velocity, respectively [37]. The Fermi energy and the carrier mobility of the monolayer graphene are assumed to be ${E_{\rm{F}}}{\rm{=}}0.6{\rm{eV}}$ and $\mu=1.3793\times{10^{4}}{\rm{c}}{{\rm{m}}^{\rm{2}}}{{\rm{V}}^{{\rm{-1}}}}{{\rm{s}}^{{\rm{-1}}}}$ [38, 39], and correspondingly the relaxation time is $\tau=0.8276{\rm{ps}}$. The refractive index of the silicon is ${n_{{\rm{Si}}}}=3.4164$. The corresponding dispersion diagram of the SPP on the air-graphene-Si interface is shown in Fig. 1(b).

Here we consider that the incident wave is at 5.368THz and the corresponding wavelength is $\lambda=55.9{\rm{\upmu m}}$. The incident angle of paraxial beam is $\theta={44.427^{\circ}}$. The thickness of grating is $d=5{\rm{\upmu m}}$ and the period of grating is $\Lambda=4.101{\rm{\upmu m}}$ for the phase matching condition. We numerically calculate the spatial spectral transfer functions ${H_{xx}}$ and ${H_{xy}}$ for the metasurface. The requirement of $\left|{{A_{x}}}\right|=\left|{{A_{y}}}\right|$ can be achieved by appropriately designing the geometric parameters of the gratings, with the width of grating slot $W=1.1323{\rm{\upmu m}}$ and the height of grating slot $H=0.8194{\rm{\upmu m}}$.

Fig. 2(a) and (b) show that the distributions of $\left|{{H_{xx}}}\right|$ and $\left|{{H_{xy}}}\right|$ are mainly dependent on ${k_{x}}$ and ${k_{y}}$, respectively. Furthermore, the blue (green) line in Fig. 2(c) and (d) corresponds to the amplitude and phase of ${H_{xx}}$ (${H_{xy}}$) along ${k_{y}}=0$ (${k_{x}}=0$). As expected by Eq. (1) and Eq. (2), Fig. 2(c) shows that $\left|{{H_{xx}}}\right|$ and $\left|{{H_{xy}}}\right|$ exhibit good linear dependency on ${k_{x}}$ and ${k_{y}}$, and the minimums occur at ${k_{x}}=0$ and ${k_{y}}=0$, respectively. More importantly, they almost overlay with each other, which indicates that the requirement of $\left|{{A_{x}}}\right|=\left|{{A_{y}}}\right|$ is approximately achieved. Fig. 2(d) shows the $\pi$ phase shift at the point ${k_{x\left(y\right)}}=0$ for ${H_{xx}}$ and ${H_{xy}}$, which is required by the first-order differentiation operation. In addition, the invariant phase difference between ${H_{xx}}$ and ${H_{xy}}$ is $\Delta\phi=\arg\left({{A_{x}}}\right)-\arg\left({{A_{y}}}\right)=5.4$, which can be compensated by the wave plate with appropriate thickness .

To illustrate the divergence computation, we consider a series of cylindrical vector (CV) fields, which have the same intensity distribution but different distributions of polarization states as shown in Fig. 3(a-c). Such cylindrical vector fields can be generated by an interference method [40], with two Laguerre-Gaussian beams with opposite topological charges and opposite handedness of circularly polarization:

$${{\bf{E}}_{{\bf{i}}}}=\frac{1}{{\sqrt{2}}}\left({{\rm{LG}}_{0}^{-n}{e^{-i{\varphi_{0}}}}{{\bf{e}}_{\bf{+}}}+{\rm{LG}}_{0}^{n}{e^{i{\varphi_{0}}}}{{\bf{e}}_{-}}}\right)$$

(5)

Here ${\rm{LG}}_{0}^{n}$ corresponds to the Laguerre-Gaussian beam with the electric field distribution ${\rm{LG}}_{0}^{n}={\left({\sqrt{\frac{{2e}}{{\left|n\right|}}}\frac{r}{{{w_{0}}}}}\right)^{\left|n\right|}}{e^{in\varphi}}{e^{-{{{r^{2}}}\mathord{\left/{\vphantom{{{r^{2}}}{w_{0}^{2}}}}\right.\kern-1.2pt}{w_{0}^{2}}}}}$, $r$ and $\varphi$ are the polar coordinates of the cross section of beam. The ${\varphi_{0}}$ is an initialized phase of the Laguerre-Gaussian beam. The unit vectors ${{\bf{e}}_{\bf{+}}}$ and ${{\bf{e}}_{\bf{-}}}$ are the left- and right-hand circularly polarized states of the two Laguerre-Gaussian beams respectively. When the topological charge is fixed, the cylindrical vector fields of Eq. (5) have the same amplitude distribution $\left|{{{\bf{E}}_{{\bf{i}}}}}\right|={\left({\sqrt{\frac{{2e}}{{\left|n\right|}}}\frac{r}{{{w_{0}}}}}\right)^{\left|n\right|}}{e^{-{{{r^{2}}}{\left/{\vphantom{{{r^{2}}}{w_{0}^{2}}}}\right.}{w_{0}^{2}}}}}$. Here we first choose the topological charge in Laguerre-Gaussian beam as $n=1$, and the radius of the beam waist is ${w_{0}}=3.82{\lambda}$. The amplitude distribution of incident field $\left|{{{\bf{E}}_{{\bf{i}}}}}\right|$ is shown in Fig. 3(a-c). On the other hand, the polarization state of the CV fields is only dependent on the azimuthal angle $\varphi$:

$${\bf{e}}\left(\varphi\right)=\cos\left[{\left({n-1}\right)\varphi+{\varphi_{0}}}\right]{{\bf{e}}_{\bf{r}}}+\sin\left[{\left({n-1}\right)\varphi+{\varphi_{0}}}\right]{{\bf{e}}_{{\upvarphi}}}$$

(6)

where ${{\bf{e}}_{\bf{r}}}$ and ${{\bf{e}}_{{\upvarphi}}}$ are the orthogonal bases associated with the polar coordinates $\left({r,\varphi}\right)$. Therefore, in the case of $n=1$, the polarization state ${\bf{e}}$ has a uniform angle ${\varphi_{0}}$ from the basis ${{\bf{e}}_{\bf{r}}}$, shown as the arrows in Fig. 3(a-c), with different initialized phase factors ${\varphi_{0}}={0^{\circ}}$, ${45^{\circ}}$, and ${90^{\circ}}$, respectively. We note that at the center of the beam, the orientation of resultant vector field is undefined and the amplitude of the electric field vanishes, which corresponds to a vector point singularity (V-point) [41, 42, 43].

We numerically simulate the output beams after analyzing the $x$ polarization component of reflected light by the metasurface. Fig. 3(d-f) correspond to the simulation results of the output field for the incident CV fields in Fig. 3(a-c). According to Eq. (5), the ideal result of the divergence computation is $\frac{{2\sqrt{2e}{A_{x}}}}{{{w_{0}}}}\cos{\varphi_{0}}\left({1-{{\left({{r{\left/{\vphantom{r{{w_{0}}}}}\right.}{{w_{0}}}}}\right)}^{2}}}\right){e^{-{{{r^{2}}}{\left/{\vphantom{{{r^{2}}}{w_{0}^{2}}}}\right.}{w_{0}^{2}}}}}$, which has a zero value at $r={w_{0}}$, and agrees well with the simulation results shown in Fig. 3(d-f). Although the intensities of the incident vector fields are the same, the divergence computation results exhibit different amplitudes, dependent on the polarization distribution of the incident field. Fig. 3(g) shows the numerical simulation result (red circles) and the ideal one of the divergence computation (blue solid line) around the beam centers, which agree well with high accuracy. Especially in this case of ${\varphi_{0}}={\rm{90}}^{\circ}$ [Fig. 3(c)], the divergence is absent since the incident vector fields are purely rotational.

Conventionally, the concepts of source and sink are used to describe the polarization distributions in vector field, respectively, with positive and negative discrete charges. By computing the divergence, the discrete charges should redistribute themselves in a continuous manner. Indeed, the intensity image obtained by our device, shown in Fig.  3(d-f), present the density distribution of the equivalent charges in vectorial paraxial beams.

The density of equivalent charge distribution exhibits more complicated forms if topological charge of incident Laguerre-Gaussian beam has higher orders. Fig. 4(a)-(c) correspond to the cylindrical vector fields of Eq. (5), with different topological charges $n=2$, $3$, and $4$. The arrows in Fig. 4(a-c) exhibit the distributions of the polarization states. The sources and sinks of the polarization distributions are represented by the sign $\pm$ in Fig. 4(a-c). Fig. 4(d-f) are the simulation results of the optical computation of the divergence, performed by the proposed metasurface. In a good agreement with the polarization states as shown in Fig. 4(a-c), the amplitude of the output field quantitatively models the density distribution of positive and negative charges. More interestingly, Fig. 4(d-f) intuitively show that out of central area (the dashed circles) there are opposite charge distributions, which is not so obvious to infer from the polarization states in Fig. 4(a-c) directly.

We propose a periodic plasmonic metasurface to realize the optical computation of divergence operation. Moreover, we demonstrate the divergence operations for a series of cylindrical vector fields, which extract the equivalent charge distributions with different polarization states or topological charges. Our study shows that such optical computing devices provide a new direct pathway to elucidate specific polarization singularities in vector fields. It would be more helpful in the three-dimensional cases, which require much more stable measurement [44, 45]. Also such real-time high-throughput analyzation of vector field provides possible applications in microscopic imaging, target recognition, and augmented reality.

The authors acknowledge funding through National Natural Science Foundation of China (NSFC Grants No. 91850108 and No. 61675179), National Key Research and Development Program of China (Grant No.2017YFA0205700), the Open Foundation of the State Key Laboratory of Modern Optical Instrumentation, and the Open Research program of Key Laboratory of 3D Micro/Nano Fabrication and Characterization of Zhejiang Province.