Inverse-designed gyrotropic scatterers for non-reciprocal analog computing

2.1 Ensemble of gyrotropic particles

The geometry for the proposed analog computing framework is described by a two-dimensional problem in which we assume that the fields are transverse electric (TE), which can be reduced to a scalar Helmholtz equation in terms of the magnetic field component H _z . We compose the scatterer by a system of small particles, where each particle is approximated by a dipole moment.

The polarizability of each particle is derived from the dominant scattering coefficients. We consider a single circular scatterer in a vacuum with radius a. The solution to an incident plane wave with e ^iωt dependency reads

where J _n and H n ( 2 ) are the nth order Bessel function of the first kind and Hankel function of the second kind, respectively, k ₀ is the free space wavenumber, θ is the angle of the incident planewave, and r and ϕ are the polar coordinates [23]. The scattering coefficients S _n are obtained by applying the interface condition at r = a which is described by

where n ̂ is the interface normal and the L, R subscripts indicate the left and right elements of the interface. For perfect electric conductors, they are [23]

Scattering coefficients for both gyrotropic and dielectric media can be described by the same expression [23]

In the above, ɛ ₀ and ε e ff = ε 0 ε t 2 + ε g 2 ε t are the free-space and effective permittivity, respectively, where ε g = ω p 2 ω c ω ω c 2 − ω 2 and ε t = 1 + ω p 2 ω c 2 − ω 2 are components of the permittivity tensor for a gyrotropic medium (Appendix A.1). Here, ω the field angular frequency, ω _p and ω _c the angular electron plasma and cyclotron frequencies, respectively [24]. Scattering coefficients for the reciprocal scenario are recovered when ω _p = 0.

The scattering coefficients S _n of order −1, 0, and 1 dominate the responses when the radius is much smaller than the background wavelength, i.e., a ≪ λ ₀. Then, we can approximate the scatterer by dipole moments, where the relation between their amplitude d = [m, p] ^T and the incident field F i n c = [ H i n c , E i n c ] T is expressed as d = αF ^inc, or

Here, m and p are magnetic and electric dipoles, respectively, and α is the polarizability tensor. Without loss of generality, due to the small size of the scatterers, no coupling between H and E is assumed, i.e., α ⁰¹ = α ¹⁰ = 0. Then, as shown in Appendix A.2, we have

In the above, we have S ₋₁ = −S ₁ for a dielectric case, which implies reciprocity.

According to [25], the magnetic and electric field induced by a dipole can be expressed by F ^scat = Γd or

where Γ is the dipole-dipole interaction matrix that maps a dipole to a scattered field (Appendix A.3).

The multi-scatterer generalization of eq. (5) for ith particle can be written as

Here, the effective incident wave for the ith particle F ̃ i i n c is composed of the incident field F i i n c and scattered field induced by other particles Γ _ij d _j , where Γ _ij is the dipole–dipole interaction matrix in eq. (7). α _i is the polarizability tensor for ith particle and the corresponding dipole moment is denoted by d _i . Following the notation of a two-scatterer case in [26], (8) can be written as

Here, each element of D is the dipole moment d = [m p] ^T induced by each particle. Each component of the vector P corresponds to the incident field F ^inc sampled at the position of each particle:

where Ω is the domain that contains all particles and δ(x − x _i ) is a Dirac delta.

Similarly, the generalization of (7) for an arbitrary position x with an N-particle system is

where each interaction matrix Γ x x i relates the field at x due to the dipole at x _i . Rearranging and substituting (9) into D of (11), then replacing P with (10), the relationship between the scattered field to the incoming field can be obtained as

The incident field can be expanded in a functional basis of our choice [5], [12], [25], that is, expressed as a summation of basis functions with distinct coefficients. These coefficients serve as control variables encoding the input signal [5]. Then, the complexity of the problem is governed by the number of input and output ports. Let B _n (x) = J _n (k ₀ r)e ^in(ϕ−θ) represents the basis for the incident magnetic field H ^inc, or the input ports. Then, the expansion reads

where g n ( x ) = B n ( x ) , ( i ω ε ) − 1 curl B n ( x ) T . The coefficient c n i n determines the amplitude of each port, and can be varied to represent the desired input function.

Likewise, the choice of output ports is motivated by the general solution of the scattered field in an exterior problem, i.e.,

where h m ( x ) = H m ( 2 ) ( k 0 r ) e i m ( ϕ − θ ) . Let h ^m (x) denotes the dual function of h _m (x). Then, the output coefficients c m o u t of (14) are obtained by

It is noted that the choice of basis functions is, in general, arbitrary; the above selection ensures completeness for the two-dimensional problem [27].

Let Q = 1 , 0 T denotes a linear operator that extracts magnetic field H from the field vector F such that H ^scat = QF ^scat. Then, eq. (12) can be expressed in terms of input and output ports as the scattering matrix equation:

or in tensor form,

The scattering matrix S _mn is a linear operator mapping the input vector c ⁱⁿ to the output vector c ^out .

2.2 Mathematical problem

The target linear mathematical problem for any ODE can be stated as

where A is the operator acting on an unknown function u and f is a given function. The operator A is not necessarily self-adjoint, which may describe a nonreciprocal physical problem. For example, we consider the following operator, which represents a typical form of a second-order boundary value problem:

with two complex-valued parameters α and β and the periodic boundary condition over 2π. In the weak form, the above strong form can be expressed as

In the above, v denotes a test function, and V is a set of admissible functions.

Among many choices of finite-dimensional approximations, we take Rayleigh-Ritz type discretization, which reads

The expansion is evaluated at the interface with a radius of choice r = a. The above basis choice implies that we use sinusoidal basis for both input and output ports in (13) and (15). Thus, f is now expressed with a finite series of sinusoidal functions representing the input.

Using this result, eq. (18) can be approximated to a matrix equation

where u and f are the vectors of the coefficients u _n , f _n defined in eq. (21), respectively. The operator A is the finite approximation of the operator A ⋅ , whose elements are calculated using

The aim is to design a scatterer such that its geometry maps the response of eq. (22) to its scattering matrix in a way that the solution is imprinted on the scattering fields when illuminated by incident fields f. Therefore, the scattering matrix S should satisfy

where I is the identity matrix and γ is a scaling factor [5].

2.3 Scatterer design

The position of the particles can be formed as an inverse design problem. From eq. (24), it follows that the cost function for this inverse design problem can be formulated as

where ‖ … ‖ 2 2 is the square of the Euclidean norm and p is the position vector of the particles. The scaling factor is optimized in addition to the positions to introduce an extra degree of freedom in the optimization problem. A new vector can be formed q = (γ, p), which is the optimization variable. A constant number of D _c dielectric and M _c gyrotropic cylinders were chosen at the start of each optimization, the optimization variable is defined

where p j i is the position of the ith particle in the j axis and R is a positional constraint. This constraint limits the optimization region to a square region of size R × R to ensure that: (a) the scatterers remain within the simulation limits and, (b) the desired compactness is achieved. Since the multiplier is an unbounded variable, the minimum constraint on γ ⁻¹ was enforced to avoid a singularity in eq. (24), while the maximum value was empirically chosen.

A two-step optimization process is used as shown in Figure 2 and the code can be found in [28]. The framework was implemented using the open source library NLopt in Python [29]. First, an optimization is performed using the global evolutionary algorithm (ESCH) [30] followed by the gradient based method of moving asymptotes (MMA) [31]. As a global optimizer is better suited to exploring a larger part of the parameter space, introducing it at the beginning allows for better exploration of a wider region of possible solutions. However, as global optimizers often suffer from slow convergence, the gradient based optimizer allows for faster convergence to the local minima of the region that the ESCH algorithm identified. The hybrid method explores the solution space more efficiently and accurately than deploying either of the algorithms on an individual basis [32].

Figure 2:

Optimization algorithm for the positions of the cylinders.

The calculation of gradients δq, required by MMA, was performed using finite differences. Given the small number of elements in vector q, the problem did not necessitate the use of the adjoint variable method. To reduce the computational time required, the process was parallelized using mpi4py [33] to calculate each element of δq on a different simultaneous process.

Random noise was integrated in the optimization process to encourage tolerance in the solution, yielding more robust devices. The gradient based optimizer was employed 10 times with a maximum number of iterations at 1,000 per optimization session. Each time the position vector with the lowest cost function was used with added noise i.e.

where q _best is the best performing geometry from the previous iterations and Y is a random variable vector and Y _k is the kth element of Y. This can also be considered a fabrication robustness measure, with geometries who are less susceptible to this perturbation being favored as the starting point for the next optimization.