Parallel integro-differential equation solving via multi-channel reciprocal bianisotropic metasurface augmented by normal susceptibilities

Analog optical signal processing has dramatically transcended the speed and energy limitations accompanied with its digital microelectronic counterparts. Motivated by recent metasurface’s evolution, the angular scattering diversity of a reciprocal passive bianisotropic metasurface with normal polarization is utilized in this paper to design a multi-channel meta-computing surface, performing multiple advanced mathematical operations on input fields coming from different directions, simultaneously. Here, the employed ultra-thin bianisotropic metasurface computer is theoretically characterized based on generalized sheet transition conditions and susceptibility tensors. The operators of choice are deliberately dedicated to asymmetric integro-differential equations and image processing functions, like edge detection and blurring. To clarify the concept, we present several illustrative simulations whereby diverse wave-based mathematical functionalities have been simultaneously implemented without any additional Fourier lenses. The performance of the designed metasurface overcomes the nettlesome restrictions imposed by the previous analog computing proposals such as bulky profiles, asserting only single mathematical operation, and most importantly, supporting only the even-symmetric operations for normal incidences. Besides, the realization possibility of the proposed metasurface computer is conceptually investigated via picturing the angular scattering behavior of several candidate meta-atoms. This work opens a new route for designing ultra-thin devices executing parallel and accelerated optical signal/image processing.


Introduction
The study of analog computing scheme has been one of the hottest areas of the optical domain over the past few years and despite the prevalence of traditional digital integrated electronic-based technologies, it has had a profound impact across the modern photonic requisitions [1]. Optical signal processing brings together various fields of optics and signal processing to acquire high-speed mathematical functions that can potentially operate at the line rate of fiber optic communications [2]. These mathematical operators can be either simple like the first-order integrator/differentiator or complex like integro-differentiation equation solvers and image processing functions.
Constant-coefficient ordinary differential and integro-differential equations represent a vast variety of basic engineering systems and physical phenomena in virtually any field of science and engineering, including temperature diffusion processes, physical problems of motion subject to acceleration inputs and frictional forces, time-domain circuit responses and so on [3]. In order to solve differential equations in photonic domain, exploiting from an all-optical differentiator or integrator is indispensable. The recent efforts to realize integrodifferential equation solver operators unveil three contemplative solutions in photonic territory. The first one requires an optical feedback loop [4,5], while the second one relies on an engineered temporal impulse response [6]. Fiber gratings [5], silicon micro-ring resonator [7], and all-optical differentiator [6] have been proposed as some candidates for implementation of these solutions. However, in addition to involving either a redundant loop [8] or an additional optical pump [9], the configuration of these schemes suffer from microelectronic limitations regarding operational speed, power consumption and significantly larger size, which is inappropriate in the new generation of optical systems [3,10]. Besides, the proposals in [4,5] only investigate the first-order differential equations. The newly sprouted terminology, i.e. 'metamaterial-assisted optical signal processing' or 'computational metamaterials' pioneered by Silva et al [11], has offered the third solution through providing a new venue for executing the signal processing functions in optics. The prominent wave-matter-interacting capabilities of metamaterial family suggest a variety of intriguing wave-based phenomena such as all-optical logics [12], all-optical differentiators and integrators [13][14][15] and some more sophisticated computing operators, such as differential and integro-differential equation solvers or/and image processing tools at the miniaturized hardware level [4,6,7,9]. The design of these freshly germinated meta-computing devices are still ongoing within two distinct frameworks of Fourier method and Green's function (GF) approach [16]. Unlike bulky designs arising from additional Fourier sub-blocks in the former case, in the GF strategy, the computing metamaterial directly realizes the desired spatial impulse response at the output without any additional subblocks [17][18][19]. Nevertheless, most of the metamaterial-assisted proposals for solving integro-differential equations are based on the bulky Fourier designs [3,10]. For instance, Zhang et al proposed and designed the dielectric metamaterial blocks to solve some particular forms of the ordinary differential equations [3]. In [10], by employing Fourier metalenses, a highly efficient but bulky metadevice was designed to get the solution of ordinary differential and integro-differential equations with constant coefficients. Albeit, there are metamaterial designs based on GF approach proposed for image processing [18]. A multitude of attempts have been also devoted to advance the functionality of all-optical signal processing architectures where several theoretical/ experimental demonstrations have been lately published [20,21]. For instance, the authors in [20] have demonstrated that the Laplacian operation can be implemented by photonic crystal slab device in the transmission mode. Zhu et al have shown that the interference effects associated with surface plasmon excitations at a single metaldielectric interface can perform spatial differentiation and edge detection without any Fourier lens [21]. Roberts et al have provided a framework to illustrate the capacity of spatially dispersive modes in uniform metasurfaces for modifying or enhancing the input image contrast and manipulating the spatial content of the incident field in analog optical information processing [22]. In a nutshell, however, all existing optical processing proposals still suffer from slow responses [11], supporting only the even symmetry operations for normal incidences [20,21], working for a certain incident angle or polarization [16,21], and solving only one integro-differential equation [3,10]. However, to the best of authors' knowledge, realizing such advanced optical signal processing functions via the GF approach has not been addressed, yet.
In another context, metasurfaces are known as two-dimensional versions of the metamaterial family [23] at subwavelength scales which cover a broad scope of applications due to their strong interactions with light and ease of fabrication [24]. They provide prominent abilities for polarizing [25][26][27][28], absorbing [29][30][31][32], channeling [33][34][35][36][37][38][39][40], and processing [41][42][43][44] waves, especially within optical frequencies. Several analytical frameworks have been recently developed to synthesize metasurfaces at microscopic to macroscopic scales such as impedance/ admittance matrices [45,46], generalized sheet transition conditions (GSTCs) by susceptibility tensors [47,48], polarizability tensors [49][50][51], and local complex scattering parameters [39,[52][53][54]. These analytical or at least semi-analytical frameworks have common basis in concept of modeling metasurfaces as zero-thickness bianisotropic sheets, gathering effective tensorial parameters. Seamlessly, they mostly ignore the presence of normal polarizations (or currents) with respect to the metasurface plane [55]. This is while adding normal polarization dramatically expands the scope of available wave-based functionalities and plays a vital role, as a new degree of freedom (d.o.f), in engineering the angular scattering behavior of the bianisotropic metasurfaces especially, when an asymmetric angular scattering behavior is of interest [55]. While such unprecedented d.o.f provided by normal polarizability components might be an avenue to mutate the properties of all-optical analog computers (like integro-differential equation solvers or image processors), there has been no report manifesting the unique scattering features of bianisotropic metasurfaces in wave-based analog computers, so far.
The problem in this paper is how to design 'computational metasurfaces' for the first time to provide suitable angular scattering diversity for realizing multiple advanced optical signal processing operators at the same time. We will mainly epitomize the operators associated with solving the integro-differential equations and applying the image processing. The basic problem can be formulated quite simply: In this article, using the GSTC/ susceptibility tensor formalism and the GF approach, we design a single reciprocal passive bianisotropic metasurface augmented by normal polarization to apply parallel/independent mathematical functions on multiple input signal fields coming from different directions. Here, we focus on the transmission operational mode as it allows the metasurface to be directly exploited as the first layer of an image recognition system. We will break the irritating restrictions of the current analog computing schemes such as implementing only single mathematical operation or even functions of transverse wave momentum at normal direction arising from angular symmetry of reflection/transmission coefficients around θ=0. Besides, the realization possibility of the proposed metasurface computer is conceptually investigated via picturing the angular scattering behavior of several candidate meta-atoms. In order to clarify the concept, we present several illustrative examples whereby diverse wave-based mathematical functionalities including solving asymmetric integro-differential equations, edge detection and image blurring have been demonstrated. The parallel metasurface-assisted scheme paves the way towards the implementation of analog mathematical operators to accelerate the optical signal and image processing routines, significantly.
2. Theoretical framework 2.1. Integro-differential equations Figures 1(a), (b) graphically sketches the main idea for performing the parallel optical signal processing based on the susceptibility-GSTCs formalism. Here, we exploit a reciprocal passive bianisotropic metasurface computer consisting of a periodic/uniform array of polarizable meta-atoms in z=0 plane, which acts as a twodimensional EM discontinuity upon illuminating by external sources. Throughout the paper, we will also assume the time-dependence w e j t omitted for conciseness. Before delving into the mathematical presentations, the main target of the paper is conceptually epitomized, here: the metasurface computer designed in this paper should provide different independent processing channels whereby it will be responsible for manipulating the transverse field profile of the input signals impinging from each channel and then reflecting/transmitting them into the same spatial channel. Indeed, the metasurface computer with a set of surface susceptibility components c c c c ( ) , , , can be thought of as the input and output signals of our linear system, where, the angular EM response of the metasurface computer (transmission or reflection coefficient versus incident wave angle) determines the corresponding transfer function in the spatial Fourier domain, respectively. Since solving the integro-differential equations in the field of optical computing has been oppressed in recent years, among all possible wave-based processing operators, our focus in this paper is mainly concentrated on the solutions of the constant-coefficient ordinary differential equation (CODE) and integrodifferential equation (CIDE). Such solutions are of great significance as these equations can model various fundamental physical phenomena and engineering systems including physical problems of motion subject to acceleration inputs and frictional forces, temperature diffusion processes, time-domain circuit responses and so on [3,10]. Regarding the input and output transverse field profiles, we consider the following second-order CODE and the first-order CIDE as [10] x e x e x d d in which, α, β, ξ, η, and τ stand for any positive parameters which can be arbitrarily chosen. We should remark that the input field, . Taking the Fourier transformation from equations (1), (2) manifests the required impulse response of the computing system in the Here, the solutions can be calculated by in which F and F −1 represent the operation of Fourier and inverse Fourier transform, respectively, and k x denotes the spatial frequency variable in the Fourier space. From a general viewpoint, the transfer function required to solve equations (1), (2) has an asymmetric behavior (no odd-or even-symmetry) with respect to k x variable. Accordingly, the current GF-based analog computing proposals fail to solve these equations for input fields coming from boresight direction, as they can realize transfer functions with only even-symmetry mathematical operators for normal incidences [20,21]. There are three important points: Firstly, the GF-based CODE and/or CIDE solving optical devices have not been reported, yet. Secondly, the current capabilities for realization of GF-based computational metamaterials [20,21] can only deal with particular cases of CIDE and CODE, albeit, for oblique processing channels. Thirdly, finding suitable meta-atoms to realize such complex spatial transfer functions will be very sophisticated and time-consuming if we deal with the problem by trial and error. Therefore, we have to provide a general synthesis framework whereby the meta-atoms occupying the metasurface computer can be elaborately tailored so that the angular trend of their reflection/transmission coefficient mimics the required spatial dependency of equations (3), (4) For this purpose, we borrow the irrefutable capacity of the longitudinal polarizability components describing the bianisotropic metasurfaces to engineer the angular scattering response of meta-atoms upon illuminating from different directions [55].

Susceptibility-GSTCs formalism
Theoretically speaking, we will derive rigorous expressions in this section to extract the necessary conditions through 'which the k x -dependency of the operators of choice, i.e. the spatial transfer functions of equations (3), (4) can be physically realized. Commonly, the GFs provided by an array of polarizable meta-atoms has a tensorial format for TE and TM polarizations which can be written as TE TM  TE TE  TE TM   TM TE  TM TM where ( ) S k x refers to the functionality of the reflection (R) or transmission (T) coefficient of the bianisotropic metasurface from the incident wave angle. The first and second superscripts also represent the polarization of the input and output waves, respectively. According to equation (5), the GF relates the transverse profile of the input field to that of the reflected or transmitted field. The spatial transfer functions belonging to the metasurface computer can be extracted using the GSTCs formalism in which the metasurface transition conditions read [47] ẃ Here, ò 0 and η 0 are the permittivity and the intrinsic impedance of the freespace, respectively. In the most general case, each susceptibility tensor appearing in equations (8), (9) includes both longitudinal and tangential components, i.e. 36 susceptibilities. Considering all of them during the metasurface design causes the problem to be practically difficult to solve, especially in the presence of spatial derivatives in equations (6), (7). Consequently, the normal susceptibilities are typically not taken into account for simplicity, as justified by the fact that it is generally possible to find an equivalent metasurface with purely tangential polarizations [56]. Such an assumption, however, has its own restrictions in our proposal. Most importantly, equating the metasurfaces having purely tangential polarizations with those exhibiting both normal and tangential ones is generally limited to the case of identical single excitation [55]. It follows that if the incidence angle of the input processing channel varies, e.g. multi-channel metasurfaces in this paper, the purely tangential metasurface will not generate the same angular scattering as the equivalent metasurface with purely tangential components. Therefore, the presence of normal polarizations plays an important role throughout this paper and bring about additional functionalities, as will be shown thenceforth. To have an intuitive awareness of how the angular scattering response of a bianisotropic metasurface depends on its normal susceptibilities, we begin the study with a simplified but yet pertinent scenario. For the sake of brevity, the TMpolarized illumination is only considered here.   [47]. In this case, the reflection and transmission coefficients in the spatial Fourier space can be expressed in terms of the metasurface susceptibilities as shown in [55], and are given by  Equations (12), (13) are spatially dispersive mathematical functions that can potentially imitate the same spatial transfer functions of equations (3), (4) in the transmission or reflection mode, respectively. Although there exist several reports on the computing systems realizing operators with odd symmetry for oblique incidences [16,18], most of the current GF-based analog computing proposals can only construct transfer functions with even symmetry at normal incidences [57].
For instance, if solving an asymmetric CIDE or CODE in the transmission mode of a normally oriented processing channel is of interest, the bianisotropic meta-atoms must treat differently against input fields propagating along θ and −θ directions. As can be noticed from equation (12), arming the metasurface computer with c xz ee and c xy em components, affiliates the transmission transfer function to k x variable, as an asymmetric function of θ. To visualize the concept, let us consider the diagrammatic representation of figure 2 in which the normal susceptibilities of the metasurface makes  T 1 2 to be intrinsically different with  T 1 3 . The other normal susceptibility components may not be effective to break the angular symmetry of the transmission transfer function, e.g. c zz ee reveals k x 2 in the relations which in turn is an even-symmetric function of θ. We can conclude that existing a normal polarization, specifically c xz ee , is inevitable for breaking the compulsive angular symmetry accompanied with the transmission function of the bianisotropic metasurfaces. Eventually, the GF-based GSTCs formalism may feasibly offer a zero-thickness homogeneous, reciprocal, and passive surface for solving the general form of the integro-differential equations in the transmission channel, without any additional bulky Fourier lens.
From a realization point of view, the angular response of the individual meta-atoms should be conceptually assessed, because the angular scattering properties of a given reciprocal metasurface has the same symmetries as those of its meta-atoms [55]. Geometrically, three types of symmetries are conceivable for the constituent metaatoms [58]: a reflection symmetry through the z-axis (σ z ), a 180°-rotation symmetry around the y-axis (C 2 ), and a reflection symmetry through the x-axis (σ x ). Regarding a reciprocal metasurface, the angular spectrum of the reflection coefficient exposes a σ z symmetry [58], while that of the transmission coefficient has a C 2 symmetry [55,58]. Most importantly, a bianisotropic metasurface with both normal and tangential susceptibilities is macroscopically achieved when the occupying meta-atoms do not microscopically render any geometrical symmetry. In this case, in line with our intention, the angular scattering response does not append any symmetry beyond those dictated by reciprocity. In fact, when the meta-atoms are deprived of any geometrical symmetry, the resultant metasurface will present c xz ee , c zx ee , c xy em and c yx me were required for realizing the asymmetric mathematical operators. In order to demonstrate the concept, let us consider the meta-atom of figure 3(a) comprising of two identical gold rods located a the distance 130 nm from each other. They have a square cross section of 380 nm and a length of 1.67 nm and are periodically repeated with the periodicity of 2 μm and 1.28 μm along the x-and y-directions, respectively. The presented meta-atom has ideal symmetry through both x-and z-directions so that it reveals all C 2 , σ z and σ x symmetry properties. The meta-atom of figure 3(d) consists of a gold rod and a L-shaped structure with the interspace distance of 130 nm. It has the same square cross section and the same periodicity but with the arm length of 1.06 nm. Contrarily, such a scattering particle would neither be σ z symmetric nor C 2 symmetric but it would be σ x symmetric. The phase and amplitude of the resulting transmission (solid lines) and reflection (dashed lines) coefficients are plotted for both meta-atoms in figures 3(b), (c), (e), and (f), respectively. As can be deduced from figures 3(b), (c), the transmission amplitude and phase results of the first meta-atom exhibit a perfect angular symmetric behavior. As expected, the second meta-atom breaking both vertical and horizontal mirror symmetries, has desirably an angular asymmetric transmission (see figures 3(e), (f)). Conceptually, the proposed meta-atom of the second type can be potentially established to elaborately mimic the required k x -dependency (non-locality) of equations (3), (4)

Metasurface-assisted optical signal processing
Up to now, we have proposed the concept of a bianisotropic metasurface with normal polarizations potentially capable of realizing the GF (kernel) associated with general linear mathematical operations. The problem of designing a metasurface exhibiting desired k-dependent (θ-dependent) non-local transmission or reflection coefficient for all transverse momentums in the range of | | k x <k 0 may not admit a general closed-form analytical solution [11]. Indeed, it is generally challenging to extract the susceptibilities of a bianisotropic metasurface to produce a desired angular transmission or reflection spectrum for performing a signal processing mission. Therefore, the metasurface synthesis will be treated as a nonlinear optimization problem which should be numerically solved. According to the role of each susceptibility tensor, for each example, the neutral components are set to zero. Then, we seek for those susceptibilities for which the difference between desired GF (for instance, equations (3), (4) and the synthesized one (transmission or reflection coefficient of the metasurface in equations (12), (13)) is minimized within the range of | | k x <k 0 . The error function is defined as the sum of squares of the differences at a finite number of N samples for both real and imaginary parts: Here,H meta (k x ) refers to the angular-dependent reflection or transmission coefficient of the bianisotropic metasurface obtained from equations (12), (13) and w re and w im represent the weight coefficients that can be established to selectively fine the real and imaginary parts, respectively. For single-operator scenarios, the optimization is applied on a certain part of the angular spectrum, while, for multi-operator metasurfaces, the best susceptibilities should be searched among all possible solutions to minimize the error for two or more distinct angular regions. Besides, it should be noted that while there are numerous combinations of susceptibility components which can be the possible solutions of the same desired Greenʼs function, most of them turn out to be unpractical. We are thus generally left with only a few practically valid solutions.

Single-channel metasurface
In this section, we will demonstrate some possible wave-based functionalities that can be unlocked by the proposed metasurface computer, as the fundamental block realizing the operator(s) of choice. One of the most important mathematical functions that has been less explored in the literature is the first-order differentiation operator (d/dx) whose spatial transfer function can be written as = ( ) G k jk x x . The reason behind this carelessness is mostly related to the complex implementation of such an k x -odd-symmetry transfer function [19], especially for normally incident input signals. It should be remembered that, in this case, an asymmetric transmission coefficient with respect to the incident angle must be presented by the designed metasurface computer. The full-tensorial description of the analog computing interface in this paper gives us an unprecedented d.o.f to masterly design a simple reciprocal bianisotropic metasurface with longitudinal susceptibility components applying d/dx operator on both polarized and un-polarized normal incidences. Figure 4 displays the synthesized and required GF of the designed metasurface programmed for transmitting the first-order derivation of the input field-profile.
As can be seen, the designed interface having normal polarization successfully mimics the exact asymmetric spatial transfer function with an odd phase jump at θ=0. The corresponding susceptibility components are given in the caption of the same figure. To inspect the performance of the metasurface differentiation, it is normally excited by different input spatial signals (see figure 5(a)) and the transmitted signals are presented in figure 5(b). In an exact accordance with our expectation from the first-order differentiation, the cosine function is converted to the negative sinusoidal profile, the Gaussian-shaped signal is mapped to its exact derivative, the vertical jumps in the input signal cause delta-like signatures in the transmitted profiles, and finally spatially constant regions in the input field yield a zero zone in the output signal. In overall, the exact 1st-derivative of all input signals has a perfect overlap with the simulated one egressing the metasurface from normal transmissive channel. This is a basic mathematical operator which will be expanded in the next section for solving CODEs or CIDEs.

Multi-channel metasurface
Smart realization and parallel computing can be pointed as two general and effective methods to accelerate the signal processing schemes. The parallel processing is aimed at duplicating the computational channels of a single metasurface to execute multiple different mathematical operations at the same time. Each processing task is performed entirely and separately by a certain spatial channel opened by the metasurface computer. Merging all required functionalities within a metasurface-based architecture, in turn, will reduce the demand for additional hardwares, enhance the analog computing rate, and the compatibility with integrated photonic meta-devices [3,19]. The hardware acceleration is one of the key necessities for todays modern optical processing, especially in commercial applications with the need for real-time performance. In this section, we explore a multi-operator bianisotropic metasurface, which enables independent parallel channels for solving the integro-differential equations, as the focal spot of this paper.
Consider a computational scenario in which two input signals collide to a reciprocal bianisotropic metasurface from two input processing channels with distinct angular directions. Upon disparate interacting with the metasurface computer and depending on the k x -modulation of the spatial transfer function dictated by the surface susceptibility tensor, the input signals of different channels may be differently mapped to the corresponding transmitted field profiles. Conceptually, the schematic sketch of the multi-channel analog computing metasurface is shown in figure 6. To demonstrate the parallel computing capability of the proposed scheme, we design a bianisotropic metasurface integrating two separate but simultaneous processing channels, one of which is dedicated to solve the k x -asymmetric second-order CODE of equation (1), for normal input signals (θ i =0°) in transmission mode and the other one is devoted to the 1st-differentiation operation for the oblique input signals (q ¹  0 i ) in reflection mode. As discussed,˜( ) G k x 1 associated with the CODE solving operator has generally an asymmetric behavior around k x =0. Therefore, we must exploit normal susceptibility components to implement such an operator. The CODE's channel parameters are set as θ i =50°, β=1.3 and α=1. For the sake of comparison, figures 7(a), (b) plot the exact and the synthesized spatial transfer functions (the phase and amplitude of the transmission or reflection coefficient versus the incident angle) associated with these two channels opened by the designed metasurface. Figure 7(a) is related to the second-order CODE in transmission while figure 7(b) pertains to the 1st-differentiation operation in reflection. The constituent susceptibility components involved in the metasurface design are listed in caption of the same figure. As can be noticed, relying on the unique features arising from those normal susceptibility elements required for imposing the k x -asymmetric behavior, the transfer functions associated with the both operators of choice are effectively synthesized. Arbitrarily, the input signal of both normally and obliquely-oriented channels are chosen as a Sincshaped function. Figures 7(c), (d) present the input and output signals of the normally and obliquely oriented channels, respectively, as well as the exact solutions for the considered CODE and the 1st-differentiation. We can clearly see that the transmitted and reflected signals captured at θ i =0°and θ i =50°directions comply well  (2) has an asymmetric behavior around k x =0. A multi-channel analogous metasurface is adopted here to solve a k x -asymmetric CIDE (θ i =0°) in the normal transmissive channel and reflect 1st-differentiation operation of the input signal coming from oblique channels (q ¹  0 i ), at the same time. In this case, we set the CIDE's channel parameters as θ i =63°, τ=1.5 and ξ=1.3. Figures 8(a), (b) plot the exact and the synthesized transfer functions in the spatial Fourier domain, describing the signal processing transformation dedicated to each computing channel. Figure 8(a) demonstrates the GF of CIDE solving operator while figure 8(b) belongs to that of the 1st-differentiation operator. Moreover, figures 8(c), (d) display the input and output signals when passing through the transmissive (θ=0°) and reflective (θ=63°) channels of the elaborately-designed metasurface. For comparison, the exact solutions of the CIDE and the 1st-differentiation are also depicted in these figures. The conformity between the simulated and exact results remains perfect as long as the spectral beamwidth of the input signals lies inside the accurate region. In conclusion, in spite of dealing with single and particular form of CODE/CIDE [3,10] in bulk scales, the zero-thickness multi-channel metasurface of figure 1 can correctly and robustly do the mission of solving multiple general (even asymmetric) CODEs/CIDEs without any Fourier lens. Eventually, it is worth mentioning that the proposed GF-based approach can be also utilized to solve higher-order differential equations with a more complex functionality from k x (supplementary appendix B).

Metasurface-assisted image processing
One of the other advanced applications of analog optical computing is image processing in which the spatial filtering is known as a fundamental step. To recognize objects inside a pre-specified image, the spatial differentiation enables the extraction of boundaries between two regions with different texture characteristics, referred to the well-known edge detection terminology [17,21,57]. On the other hand, Gaussian smoothing is an image processing function which is widely applied to the images, typically to reduce noise and detail [59][60][61]. In image blurring, the edge content is diminished to some extent, making the transition form one region to the other one very smoothly. It can be conventionally established in the vision, photography, and medical imaging applications [59]. Mathematically, applying a Gaussian blur to an image is the same as convolving the image with a Gaussian function. From the Fourier viewpoint, this is also known as a two-dimensional Weierstrass transform [62]. Keeping these two functionalities in mind, we wish to assign other processing missions to the multioperator bianisotropic metasurface: transmitting the edge-detected image and reflecting the blurred image at the same time for input images coming from normal processing channel.
The GF associated with 1D blurring process has the following spatial content: a relatively high intensity near the center (k x =0) with a gradual reduction toward the periphery whose treatment is often approximated by a Gaussian curve [60,61]. Meanwhile, the GF describing 1D edge detection processing closely resembles the 1st-differentiation operator. We now demonstrate various aspects of the parallel image processing unlocked through the proposed reciprocal passive bianisotropic metasurface having normal susceptibility components. In this way, the homogeneous susceptibilities are tailored in a manner that the angular spectrum of the reflection response belonging to the metasurface obeys the transfer function of the edge detection operator while that of the transmission coefficient approximates the Gaussian-shaped transfer function pertaining to the image blurring operator, in the vicinity of θ=0°(see figure 9(a)). Figure 9(b) compares the exact and synthesized spatial transfer functions at reflection and transmission modes. The results illustrate a perfect consistency inside the region of  | | W k 0.35 0 . The quantities of the utilized susceptibilities are quite feasible and listed in the caption of the same figure. To inquire the performance of our multi-operator metasurface in realizing parallel image processing, the image illustrated in figure 10(a) is considered as the input image field where its inner and outer domains have different intensities. The transmitted and reflected images are numerically simulated and the corresponding transverse field profiles are displayed in figures 9(b), (d), respectively. As expected, the 1D edge-detector metasurface successfully exposes all outlines of the normally incident image along the horizontal orientations, in the transmission mode (see figure 9(b)). This is while all the prevalent edge detection schemes compulsorily operate at oblique incidences [21] or need bulky Fourier lenses to do their mission upon illuminating by normally input signals [63]. Figure 9(c) demonstrates the feasibility of image blurring process, where the high-frequency components of the input image are quite removed when it is spatially mapped via the bianisotropic metasurface into the reflected image. As emphasized, the proposed bianisotropic metasurfaces can realize 1D edge detection and since the differentiation has been executed along x direction, only the vertical edges have been extracted. In a similar setup, a bianisotropic metasurface carrying out differentiation along y direction could be employed to resolve the horizontal edges. Besides, in order to extend these functionalities to 2D scenarios, two different ways may be feasible: combining two identical differentiators rotated by 90 • with respect to each other [19] and cascading two separate differentiators where the first layer is responsible for extracting edges along the x (or y) direction and the second layer is designed for detecting edges along the y (or x) direction, may be feasible [18]. In this way, all the horizontal and vertical edges of the original image can be elaborately resolved.
In overall, armed with a homogeneous reciprocal and passive bianisotropic metasurface exposing normal susceptibilities, the proposed multi-channel architecture executes the multiple delivered image processing missions, simultaneously. We again emphasize one the vital role of normal polarization which cannot be replaced with a dual tangential susceptibility profile, because the metasurface computer of this paper deals with multiple projections from different angular directions [55].

Methods
The results of figure 3 were simulated by using Time domain solver in CST full-wave commercial software. The open add space boundary conditions were assigned to the walls along all directions and a plane wave with  variable angle was utilized as the excitation source. The amplitude and phase of the scattered fields at the far-field domain have been recorded and the corresponding results have been given in figures 3(b), (c), (e), and (f). The spatially-dispersive transfer functions (GFs) were obtained by the closed-form non-local reflection or transmission coefficient of the designed bianisotropic metasurfaces in equations (12), (13) Moreover, using the Fourier expansion of the input fields as the superposition of plane waves with different illumination angles, the output fields were calculated by ò . Here, f (x) denotes the transverse profile of the input field, W is its spatial bandwidth, and k 0x = k 0 sinθ inc refers to the wavenumber of the zero-harmonic plane wave along the x direction.

Conclusion
To conclude, we have theoretically illustrate how the normal susceptibilities of a homogeneous passive and reciprocal bianisotropic metasurface can be tailored to realize different image/signal processing missions at the same time. Various appealing wave-based functionalities have been unlocked for the first time via the proposed GSTC/GF-based design framework including, but not limited to, solving CODEs/CIDEs only by a single zerothickness interface and multi-channel image/signal processing implementing asymmetric/symmetric spatial transfer functions upon normal incidences. The sophisticated realization of the latter case makes the previous works to address only the even-symmetric 2nd-differentiation in their demonstrations [20], leaving the room free for more advanced signal processing operators. The strategy presented in this paper, however, would be very promising as the basic design framework for many all-optical parallel computing architectures. Several illustrative simulations for demonstrating signal and image processing functions have supported this study and their realization possibility has been conceptually assessed. The proposed design framework overcomes the substantial restrictions imposed by previous works such as large architectures arising from the need of additional sub-blocks, slow responses, executing only single processing mission, and most importantly, supporting only the even-symmetric mathematical operations for normal input signals. We expect this work to pave the way towards realizing compact all-optical computing devices for executing multi-channel advanced signal/image processing [17,64].