Curved boundary integral method for electromagnetic fields

: The angular spectrum method is a rigorous method to synthesize near and far-field electromagnetic beams from planar field distributions. However, this limitation of planar surfaces has restricted its applicability to beams with simple focal planes. We propose a curved boundary integral method (CBIM) to synthesize electromagnetic beams from arbitrary surfaces to address this limitation and expand the method’s scope to synthesize beams from and between shaped objects. This study presents a detailed theoretical framework behind the CBIM and validates its effectiveness and accuracy with a comprehensive set of simulations. Additionally, we present mathematical proof to support our proposal. The proposed method satisfies Maxwell’s equations and significantly benefits optical systems and inverse beam design. It allows for analyzing electromagnetic forward/backward propagation between optical elements using a single method. It is also valuable for optical force beam design and analysis.


Introduction
The angular spectrum method is a rigorous method to synthesize near and far-field electromagnetic (EM) and acoustic fields based on the principles of diffraction theory.Its underlying concept stems from the fact that any of these fields can be presented as a superposition of individual plane waves, each with its unique amplitude, phase, and direction [1][2][3][4][5][6].The ASM provides a robust framework for examining wave propagation in various scenarios by decomposing the fields into constituent plane waves.
The ASM has two main advantages in electromagnetic simulation and analysis [7][8][9][10][11][12][13][14][15][16].First, due to its source-free eigenfunction nature, the incident field can be synthesized from known planar electric field distribution without singularities.Secondly, it addresses complex wave phenomena, including diffraction, interference, and scattering, often encountered in practical applications [17].However, ASM is limited to synthesizing EM fields from planar field distributions.The difficulty with handling curved, non-planar surfaces hampers its applicability in various fields that require precise control and manipulation of electromagnetic beams, such as optics, imaging, and beam-shaping applications.
This article proposes a curved boundary integral method (CBIM), which is closely related to the 2D ASM, except that the surface can be curved in the CBIM.Both methods are locally equal on the propagation axis's positive side but have some exceptions on the negative side.Thus CBIM is strictly a source method, while 2D ASM is a source-free method, allowing the beams to propagate through a field distribution.
In [18], the planar ASM was expanded to approximate diffractive fields from 2D curved surfaces.The expansion was obtained by dividing the curved 2D surface into the stepwise subregions windowed by a Gaussian function.This approach was later expanded to the 3D surfaces in [17], where the authors used the same planar stepwise subsections with Gaussian distribution.Both methods approximate well the diffraction, reflections, and transmission from the curved single-layer boundaries when the radius of curvature (RoC) is much larger than the wavelength (e.g.>100) [19].
The CBIM derivation begins with the same approach of dividing an arbitrary surface into planar subregions, presented as equal amplitude areas without Gaussian distribution or windowing (e.g.Stepwise AS).Then the areas of these subregions are shrunk into infinitesimally small source points, which approach a Dirac delta function in the limit.This approach yields exact results (in the limit of infinitesimally small partitions or infinite frequency) of the diffractive field resulting from an arbitrary surface distribution considering surface structures comparable to or larger than the illumination wavelength, i.e., kα ≥ 4π, where k is the wavenumber and α is the radius of the sphere.Addtitionally, the advantages of the presented method are non-singularity at the source points, allowing the synthesis from any compactly supported and continuous electric field distribution on a regular surface.We have utilized these results in [20][21][22][23].
Moreover, each source point is analogous to the angular spectrum presentation of the Stratton-Chu method with added polarization components.The Stratton-Chu method obtains the electromagnetic radiation from a source point as the curl of a magnetic dipole multiplied by scalar Green's function [24].This presentation satisfies Maxwell's equations and approximates well-curved surfaces when the RoC is larger than a few wavelengths [25].

Curved boundary integral method
Section 2.1 proposes a CBIM presentation for any compactly supported and continuous electric field distribution defined on a regular surface.Otherwise the distribution can be chosen arbitrarily.Section begins by introducing a general parameterization of an arbitrary surface.Then the local base vectors are derived from parametrization to determine a Cartesian coordinate system for surface sources.Further, the electric field polarization is modified by base vector rotation.Also, a transformation matrix is presented for mapping the local source coordinates in a global coordinate system and vice versa.After coordinate mapping, an arbitrary surface parameterization is discretized into infinitesimally small surface elements.The key concept is that Riemann's surface integral combines differential surface elements [26], which can be approximated as locally planar elements.The CBIM models the radiated field from each element, and the total electromagnetic field is the sum of the modeled fields from each source point by the superposition principle [1].This approach relies on the assumption that when the areas of the elements are sufficiently small, the synthesized field from the differential elements approaches the field from the original surface.Section 2.2 details a mathematical proof of the CBIM presentation.Section 2.3 expands the CBIM presentation for the magnetic field.Section 2.4 introduces the repropagation method between surfaces.

Electric field representation
Let the global Cartesian base-vectors be (e x , e y , e z ) in space R 3 , where position vector is r = xe x + ye y + ze z and global coordinates are marked as r glob = (x, y, z).Consider a regular, compact surface Ω in R 3 , which has a continuously differentiable parametrization with parameters p and q as Let us place the origin of the local coordinate system at point o(p, q) on the surface Ω, where two base vectors f 1 and f 2 define tangential plane at the point, and the third base vector f 3 is a local surface normal vector, see Fig. 1.These vectors will be nonzero due to regularity of the surfase.When normalized, these vectors form an orthonormal base (f 1 , f 2 , f 3 ) and are defined as The final local base (e 1 , e 2 , e 3 ) is obtained from the base (f 1 , f 2 , f 3 ), when the local polarization of the electric field is chosen.Let us set where the sign is chosen according to the desired electromagnetic propagation direction.In a simple case, when polarization is known on the surface Ω, the local electric field vector e 1 on Ω can be presented with the base (f 1 , f 2 ), and the local magnetic field vector e 2 is obtained as a cross-product with local propagation direction e 3 as e 2 = e 3 × e 1 .Polarization on the surface Ω can also have externally sourced initial conditions, for example, from a polarization of an incident beam, that creates the electric field distribution on the surface Ω.Let p = p(p, q) be a vector that is perpendicular to the plane of local field polarization.Then, local direction e 1 is an intersection of the plane normal to p, and the tangential plane of surface Ω spanned by (f 1 , f 2 ).Thus e 1 , the base vector along the local electric field, and e 2 , the base vector along the local magnetic field, are obtained as The vectors (e 1 , e 2 , e 3 ) form a right-handed orthonormal base, which together with the origin o defines a local coordinate system.Let vector r be presented in this base as r = xe 1 + ȳe 2 + ze 3 .Its coordinate triple is marked briefly as r loc = (x, ȳ, z).When the vector r is expressed with both (e x , e y , e z ) and (e 1 , e 2 , e 3 ) bases, as it is known, the dependence of the coordinates on each other is determined by an orthogonal transformation matrix where the columns are the coordinate triples e 1,glob , e 2,glob and e 3,glob .Let the position P be presented in global coordinate system as r = xe x + ye y + ze z .Position P presented in the local coordinate system with a global base is Let us denote briefly r = r − o.Thus, we have presentations of the location r in both local and global coordinate systems.The dependence of the corresponding coordinates can be written using the transformation matrix Eq. ( 5) as Let us divide the segment of curved surface Ω containing the electric field distribution E 0 (p, q) into a finite partition of small separate sets ∪ c Ω c = Ω.Select the point o(p, q) on differential element Ω c , the vectors associated with it are e 1 and e 2 .Area Ω c is projected to the plane spanned by the vectors e 1 and e 2 .The resulting area at the plane is marked as Ω t .Likewise, the function E 0 is projected into that plane, limited to set Ω t and zero elsewhere.This geometry is presented in Fig. 2. The resulting projection is marked as E t 0 (source point amplitude and phase).In the local coordinate system it holds z = 0 on the piece Ω t , because this is located on the xȳ−plane.E t 0 is valid when x ≈ 0 and ȳ ≈ 0. The angular spectrum theory is locally applied on the local coordinate system to function E t 0 valid on piece Ω t .Let's observe the field it creates at point r, where the coordinates are as in Eq. ( 6).At a fixed point r holds [1] where is the wavenumber.The total field from the arbitrary surface is obtained by summing the fields from the differential sources by the superposition principle as where local coordinates (x, ȳ, z) as in Eq. ( 6) and the total field is polarized along the global e 1 unit-vector.The obtained field Eq. ( 9) also satisfies the Helmholtz equation as unit vector e 1 is locally constant.To be precise, if the limit with respect to the areas |Ω t | of the sets Ω t exists, we define As noted at the beginning of this section, this limit also works physically when considering smooth surfaces.It is even more precise when the wavelength decreases compared to the local radius of the surface.Finally, a further examination is made to obtain a closed form for clause Eq. ( 10); we present proof in the next section.
Heuristically, when |Ω t | shrinks, the function normed in volume, approaches a Dirac delta δ (0,0) in the sense of the distribution theory; 0 → E t 0 ( 0, 0)|Ω t |δ (0,0) .As known, the Fourier transform of Dirac delta is 1; thus, by small since E 1 (p, q) = E t 1 ( 0, 0) when o is fixed.Then Eq. ( 9) and Eq. ( 11) yield first the form of a Riemann sum and then an integral as a limit where dΩ = ∥ ∂o ∂p × ∂o ∂q ∥dpdq.The final form gets an e 1 −component for the incident electric field and can be written as where and x, ȳ and z are as in Eq. ( 6) and notation (r; o) defines electric field on observation point r related to the location (p, q) on the parametrizied surface.An electromagnetic beam in the plane wave spectrum representation has polarization vector components along the propagation direction.To account for polarization let's expand Every single plane wave component in the integral representation of the final incident electric field, like Eq. ( 14), has to be orthogonal to k.
As in [27] we obtain where sign(z) and |z| in the exponent follows from that the propagation created at a point o is symmetrical with respect to the tangential plane, but we define Hence, a total incident electric field is Furthermore, the polarization can be defined as and we can write briefly There are two critical observations regarding this derivation.First, E 1 (r; o)e 1 is analogous to −∇ × (MG) where the scalar Green's function (G) is replaced by its angular spectrum given by Weyl identity, and where M = −2n × e 2 is a small magnetic dipole [28].Thus, Eq. ( 13) is a particular case of the Stratton-Chu equation for the electric field in which only the magnetic sources M = −n × E surf are presented and the electric ones J = n × H surf have been removed [29].Secondly, if the surface Ω is planar, Eq. ( 17) and Eq. ( 19) return at the half-space z ≥ 0 to the original two-dimensional angular spectrum method.

Proof of curved boundary integral method
Theorem 2.1 Let the function E 0 be continuous in Ω and let r ∈ R 3 be such that (r−o)•( ∂o ∂p × ∂o ∂q ) ≠ 0 for all o ∈ Ω.The field E 1 is created by Eq. ( 13) and ( 14) at point r.
Proof.Define ϵ >0.Let us show first that when the partition Ω t is sufficiently small, in other words, the areas |Ω t | are small enough regardless of t, by replacing the function E t 0 supported by each projected plane Ω t with a constant E t 0 ( 0, 0) = E 0 (p, q), the error made in Eq. ( 8) is smaller than ϵ/4.Function E 0 is uniformly continuous on the compact segment of surface Ω, and based on this, with sufficiently small |Ω t | it holds for all x, ȳ ∈ Ω t and t, where the constant I will be defined later on.In the projected plane, the Fourier transform is valid: for all , the exponential becomes real and negative.Let I be, at first, a continuous elementary function where z = z(o) = z(p, q).Then because of Eq. on the compact set {(p, q)|p ∈ [p 1 , p 2 ], q ∈ [q 1 , q 2 ]}.Moreover, |z| ≥ zmin for all t regardless of partitioning.Finally, the constant I is defined Due the monotone of integral, I(|z(p, q)|) ≤ I on all t regardless of partitions.The upper limit for the error due E 0 (p, q) on the Eq. ( 8) is obtained as where Eq. ( 20) and Eq. ( 23) have been used.Because ||e 1 || = 1 and the triangle inequality the error in Eq. ( 9) is In the above approximation, a constant E 0 (p, q) is used on a region of Ω; as a function, it is zero on the plane outside of this region.Let us denote this function simply as E 0 (p, q).Next, consider the Fourier transform of the function E 0 (p, q).Because function ȳ is integrable, there can be found an origin centered closed disk B such that for all t because |z| ≥ zmin .For all (k x, k ȳ) ∈ R 2 it holds Let us now consider the error (err 2 ) made in integral Eq. ( 8) when F {E 0 (p, q)} is replaced by the constant E 0 (p t , q t )|Ω t |.Under the estimates Eq. ( 26) and Eq. ( 27) it follows For all (k x, k ȳ) ∈ B and (x, ȳ) ∈ Ω t , we can estimate by continuity of the function always when |Ω t | is small enough, because then x ≈ 0, ȳ ≈ 0 and hence k x x + k ȳ ȳ ≈ 0 − that is the assignment of B. Then in the Fourier transform at point (k x, k ȳ) ∈ B holds Replacing by F {E 0 (p, q)} with the constant E 0 (p, q)|Ω t | in Eq. ( 8) we make an error where the estimates Eq. ( 28) and Eq. ( 30) have been used.The error for sum Eq. ( 9) is obtained as Thus, when the partition Ω t is dense enough and r ∈ R 3 be such that (r − o) • ( ∂o ∂p × ∂o ∂q ) ≠ 0 for all o ∈ Ω, we get by Eq. ( 25) and Eq. ( 32) Finally, by continuity of the integrand, the integral in Eq. ( 13) exists, and, because of the estimate Eq. ( 33), for the points r ∈ R 3 , as in Theorem 2.1, there also exists The proof above also holds for the component ∬ Ω e 3 (o)E 0 (o)E 3 (r; o)dΩ of the field in Eq. ( 17); only a few adaptions are needed.Thus, we can also obtain E 3 as a limit of polyhedra approximation and the planar theory.

Magnetic field representation
Let's fix 0 ∈ Ω and consider a point r. which has non-zero local coordinate z.Then we can take a curl of E(r; o) = E 0 (E 1 e 1 + E 3 e 3 ) at r.As well known, a magnetic field at r is where a = i/ωµ is a constant and Now the real propagation direction of Poynting vector S = 0.5Re {E × H * } in half-space z<0 can be chosen with the use of sign(z) in the H 1 and H 2 field components.With the sign(z) in H 1 and H 2 components, electromagnetic energy propagates away from the surface Ω on both ±z half-spaces.Without the sign(z) in H 1 and H 2 fields, the electromagnetic energy propagates through the surface Ω, similarly as it propagates through the focal plane in 2D ASM.If z ≠ 0, we obtain a total magnetic field at r as Simliar to the case of the electric field, we can also obtain the magnetic field as a limit of polyhedral approximation of Ω using the same project tangential plane theory on their faces at point r, where z ≠ 0 for all o ∈ Ω.The proof of Theorem 2.1 is also valid here; only a few adaptions are needed.
It can be proven that Eq. ( 17) and (37) for the electric and magnetic fields at all point r ∈ R 3 \ Ω satisfies a light condition related to Ω.For example, on sphere S 2 , all points satisfy that condition.This situation is confirmed in the simulations presented in following sections.The required extended proof is outside the scope of this work and will be presented in a separate article.Its difficulty lies primarily in that integrals do not exist in the ordinary sense when = 0.

Repropagation between surfaces
Consider the electromagnetic field (E, H) at some point g on the repropagation surface G.At that point, the direction of the field's propagation is along Poynting vector S(g) = 0.5Re {E(g) × H(g) * }.The "differential plane" T, from which the local incident field is thought to be repropagated back, is perpendicular to the direction S(g), and g ∈ T. Let's create a local coordinate system, which origin is fixed on g and an orthonormal base (f 1 , f 2 , f 3 ), where f 3 = −S(g)∥S(g)∥ −1 and f 1 and f 2 are in the differential plane T. Let's find the direction of f 1 as a real part of projection of E into the plane T as Hence, the normalized base vector for the electric field is and the second base vector for magnetic field is The complex vector E is presented in this base as and to the corresponding projection coefficient for E along f 1 is as the E • f 2 and E • f 3 components approach zero.Now we have obtained local repropagated plane wave decomposition, which is polarized along the plane (f 1 , f 3 ).Also, the transformation matrix between the global base (e x , e y , e z ) and the local base Using Eq. (43), we obtain the local coordinates at point r ∈ R 3 as Hence, the repropagated electric field is where and The repropagated magnetic field can be obtained as a curl of F(r) as in the subsection 2.3.

Electric and magnetic fields from a source point
The total electric and magnetic fields from surface Ω are obtained as a superposition of fields from source points.As shown in derivation, when the surface area of the differential source plane approaches zero, the result Eq. ( 13) is analogous to the angular spectrum presentation of , where G is a scalar Green's function and magnetic dipole is defined as M = −2e 3 × e 1 .We expanded the E 1 (r; o) presentation to include the E 3 (r; o) component, which enables propagating polarization state and allows us to compute the magnetic field with via the curl.
Next, we present the radiation patterns E(r; o) and H(r; o) from a source point, where the magnetic dipole is positioned to the origin along the y-axis.Fields are integrated over a disk B(0 }︁ to obtain propagating and evanescence fields for near-field presentation.Also, the fields are integrated over B(0 }︁ to obtain only propagating fields for far-field presentation.See the simulation arrangement in Fig. 3.The electric field polarization is now along the xz-plane, and the field propagates to the positive z-direction.Figures 4-7 illustrate the electric and magnetic field radiation with B = (0, 10k) and B = (0, k) integration limits.As seen from the simulations in Fig. 4-7) the integration limit should be selected according to the application.When the radius of curvature, of the surface on which the distribution is defined, decreases compared to the wavelength or the fields are computed close to the surface distribution, evanescent fields should be accounted for by integrating over a sufficiently wide range in the real k-domain.

Field synthesis with an increased radius of curvature
Let's examine the error created by the curved surface synthesis by comparing the results with physical optics simulations, which utilize electric and magnetic surface currents.These simulations synthesize the field from an electric field distribution on the spherical surface, polarized along the e θ −component (standard spherical coordinates).The source distribution A 0 is selected to be equal amplitude over 22.5 • solid angle with sharp edges (tophat) for creating strong interference patterns.The fields are synthesized from a sphere with the fixed radius α and evaluated on a transverse plane fixed with the size of 12.8 × 12.8 α with 200 × 200 datapoints.This plane is positioned 3.8 α distance from the sphere's origin.Then the ratio between the wavelength and radius of the sphere size parameter kα = 2πα/λ varies from 4π to 200π.The arrangement and synthesized field with the size parameter kα = 20π is illustrated in Fig. 8.The difference between the CBIM and physical optics simulations are compared as an average of difference on each datapoint on the evaluation plane as for amplitude and phase difference, respectively, where N = 200.As seen from the simulation 9, both integration limits B = (0, k) and B = (0, 10k) give reasonable accurate results when RoC≥ 2.5λ.The benefits of wider integration disk B increase when the RoC is small.We can conclude that the CBIM accurately approximates the beam synthesis from curved surfaces.

Repropagation example
To further demonstrate utility, we synthesize electromagnetic fields from a tophat electric field distribution A 0 , which is polarized along spherical e θ base vector.The A 0 is positioned on the α radius sphere's surface and propagated onto a 10 × 10α evaluation plane at the 5α distance from

Conclusions
This article introduced the curved boundary integral method (CBIM) for electromagnetic field synthesis based on the angular spectrum method.The CBIM leverages the concept of discretizing an arbitrary surface into infinitesimally small elements and models the radiated field from each of them.Simulations suggest accurate approximation of diffractive fields from surfaces with RoC larger than 2.5 times the wavelength.The CBIM's advantages include non-singularity at source points, enabling synthesis from any compactly supported and continuous electric field distribution on a regular surface inside the computational domain.
In addition, the four essential points of CBIM are that (1) CBIM is strictly a source method radiating away at both sides from the field distribution, unlike the source-free angular spectrum method.(2) CBIM is analogous to the magnetic dipole presentation of the Stratton-Chu method with polarization perpendicular to the local propagation direction.We added a polarization component to the propagation direction, which enables us to obtain magnetic field presentation via curl.Both polarization components satisfy Maxwell's equations.(3) CBIM returns to the nominal 2D ASM presentation in the half-space of positive z when the source surface is planar.(4) Each source point radiates a field to other points, modifying the total source distribution.This filters non-physical distribution E 0 to a physical one, which always satisfies Maxwell's equations.On the other hand, the final form of the more complex source distributions should be verified by simulations.
Furthermore, the article explored the repropagation of fields between surfaces, presenting a method using a local coordinate system and transformation matrix.This approach provides a comprehensive framework for simulating field interactions across surfaces.
In conclusion, the development of the curved boundary integral method extends the capabilities of the angular spectrum method by addressing its limitations related to curved surfaces.It offers new possibilities for applications in optics, imaging, and beam-shaping scenarios, where precise control and manipulation of electromagnetic beams are essential.Further research and applications of this method will advance our ability to model and simulate complex wave phenomena.

Fig. 1 .
Fig. 1.Vector relations between source point base-vectors on an arbitrary surface and a global origin.

Fig. 3 .
Fig. 3. Source point simulation arrangement, where magnetic dipole M is positioned along y-axis.The electric field is computed on xy-and xz-planes.The magnetic field is computed in xy-and yz-planes.

Fig. 8 .
Fig. 8.The simulation arrangement and global Cartesian electric and magnetic field components of the incident field synthesized from a tophat distribution with 22.5 • solid angle.In this simulation, fields are computed on a 12.8 × 12.8 α plane, located 3.8 α from the origin of the sphere with a size parameter kα = 20π.

Fig. 9 .
Fig. 9.The amplitude and phase difference between CBIM and physical optics simulations from a curved surface with RoC of α.The size parameter kα is decreased from 200π to 4π.

Fig. 10 .
Fig. 10.Electromagnetic beam synthesized from the tophat electric field distribution positioned on the α radius sphere using 40π size parameter.The synthesized field is evaluated at the 10 × 10α transverse plane positioned at 5α distance from the sphere's origin.

Fig. 11 .
Fig. 11.Amplitudes and phases of the repropagated spherical field components on the spherical surface using kα = 20π size parameter in the simulations.

Fig. 12 .
Fig. 12. Amplitudes and phases of the repropagated spherical field components on the spherical surface using kα = 200π size parameter in the simulations.