#### Abstract

In this work, the numerical steepest descent path (NSDP) method is proposed to compute the highly oscillatory physical optics (PO) scattered fields from the concave surfaces, including both the monostatic and the bistatic cases. Quadratic variations are adopted to approximate the integrands of the PO type integral into the canonical form. Then, on involving the NSDP method, we deform the integration paths of the integrals into several NSDPs on the complex plain, through which the highly oscillatory integrands are converted to exponentially decay integrands. The RCS results of the PO scattered field are calculated and are compared with the high frequency asymptotic (HFA) method and the brute force (BF) method. The results demonstrate that the proposed NSDP method for calculating PO scattered fields from concave surfaces is frequency-independent and error-controllable. Numerical examples are provided to verify the efficiencies of the NSDP method.

#### 1. Introduction

In electromagnetics (EM) community, as the product of the wavenumber and the electrical size of the considered object are large enough, the calculation of the scattered EM field belongs to high frequency scattering problems. Computing the high frequency scattered fields from electrically large perfect conducting scatterer remains significant and challenging [1–4]. In 1913, Macdonald introduced the physical optics (PO) approximation [5], which is used as an effective method for computing EM scattering problems [4, 6–8]. By this method, local tangent plane approximation is adopted to approximate the current in lit parts, and current in shadow parts is taken as zero. The PO approximation [5, 6, 8–14] is of high efficiency when calculating scattered fields from smooth and electrically large scatterer. Mathematically, the amplitude term of PO integrand is slowly varying. However, the exponential of the phase function term is highly oscillatory as the frequency *k* increases. In this sense, the computational effort of a direct numerical method [2, 11] for calculating the PO type integral is extremely high.

Some studies assume that the amplitude function and phase function terms vary linearly. Then, the physical optics integral is analytically simplified into several line integrals [15–17]. In this situation, the closed form formulas can be derived from the flat polygon patches. However, the stationary phase point (SPP), the resonance points (RP), and the vertex point (VP) contributions of the PO integral can not be captured when the phase terms of the PO type integrand are approximated into linear terms. Therefore, the traditional high frequency asymptotic (HFA) method [6, 18–21] was proposed, approximating the phase terms of the integral into quadratic forms. The HFA method could approximate the PO integrand into several dominating terms [18], corresponding to contributions from the SPP, the RP, and the VP. With the aid of HFA method, the PO-based EM fields could be calculated with frequency-independent computational workload.

The numerical steepest descent path (NSDP) approach [22–29] is proposed as an effective way to calculate the highly oscillatory PO type integral. With the aid of the NSDP method, the highly oscillatory PO type integrand is converted into several smooth ones on the complex plane by the contour deformation technique. Importantly, the contributions from the SPP, the RP, and the VP were captured through the NSDP method successfully. In former works [26, 27], the PO type integral from hyperbolic and quadratic convex surfaces has been studied. Our recent study [29] summarizes the promising progress on calculating high frequency scattered fields with the aid of the NSDP method. However, basic theories on NSDP approach have not been completed yet. Specifically, the NSDP method needs to be extended to the realistic concave surfaces.

The contributions of this work are that, first, physical optics scattered fields from concave surfaces are considered. In this case, the phase function term of the PO integrand takes the quadratic concave behavior, which differs from works in [26, 27]. Therefore, the derivations of calculating the PO type integral are different, the NSDPs are correspondingly changed, the stationary phase points and boundary resonance points are all different, and the integral paths are also changed. Furthermore, by adopting the NSDP method, the PO type integral from concave surfaces could be calculated in frequency-independent workload and error-controllable accuracy. Compared with the HFA method, the accuracy could be improved by around 2 digits. Finally, former work in [28] only calculates the PO type integral from limited concave triangular patches. In this work, however, we extend the NSDP method to the realistic concave surfaces. Together with works in [26, 27], the NSDP method could be employed to treat PO integral with all types of phase behaviors and be applied to realistic scatterer.

This communication is organized as follows. First, the PO integral under the PO approximation is derived, and the considered concave surface is discretized into several triangular patches. The Lagrange interpolation polynomial approximation and the affine transformation are used to approximate the amplitude function and phase function of the PO type integral. Next, the double integral is transformed to some highly oscillatory line integrals after integration by part. Then, we propose the NSDP method to obtain corresponding NSDPs and convert highly oscillatory PO type integrand into several smooth ones. Finally, numerical examples of calculating PO scattered fields from concave surfaces are provided, and comparisons with the brute force method and HFA method are given.

#### 2. PO Surface Integral

##### 2.1. Physical Optics Approximation

In this work, the size of considered concave surface is denoted as , the distance between the observation point and the origin is , and wavelength is . We consider the far-zone of the scatterer with and adopt the PO approximation to the accurate electric equivalent current. Then, according to the Stratton–Chu integral equation formulation [3], for the plane incident wave, the PO scattered fields can be written as

Here, is the wavenumber outside , is the incident electric polarization wave vector, is the observation point, and is the surface point on , as shown in Figure 1(a).

**(a)**

**(b)**

Equation (1) is the expression of the PO type integral of bistatic scattered electric field. The vector function in (2) is the vector amplitude function and in (3) is the phase function.

For the monostatic scattered field, with , (1) can be expressed aswith

Therefore, under the PO approximation for both bistatic and monostatic cases takes the general form [25].

The integrand contains a slowly varying amplitude term and the exponential of the phase term. When increases, the integrand is highly oscillatory due to the term, as shown in Figures 2(a)–2(d). In this sense, the computational cost of a direct numerical method grows dramatically with a larger *k*.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 2.2. Extending the Triangular Patch into the Realistic Quadratic Surface

In previous work [28], we preliminarily calculated the PO type integral from quadratic concave surfaces. However, the method had clear limitation; that is, it could only be applied for triangular patches and was not able to calculate the realistic quadratic surfaces. On the other hand, the actual geometric modeling is generally carried out by quadratic surface subdivision. In this sense, the previous work has been flawed in engineering applications, and it is necessary to develop numerical method for the realistic quadratic surfaces.

We first project the quadratic surface into the coordinate system, which is . Then, is divided into several triangular patches. In this way, the PO surface integral (7) can be expressed as

Here, is the expression of the quadratic surface . The second-order polynomials and in (7) are the approximated amplitude and phase functions, which are obtained by the Lagrange interpolation polynomial approximation [30] on the triangular patches , respectively. Formulas of and arewith coefficients .

In addition, in order to simplify the derivation, affine transformation is used to simplify the phase function into its canonical form:

#### 3. Calculating the PO Scattered Fields With the Concave Phase Function Term

##### 3.1. The Concave Phase Function Case

From (9), it can be concluded that there are four different cases of the phase function due to different combination of phase term, which are , , , and , respectively. Calculating PO surface scattered field with phase function of the first three cases has been studied by our previous works [26, 27]. In this work, we study the concave surfaces case, and a canonical PO type integral can be derived from (9):

With different phase function term of the PO integrand, derivations of calculating the integral shall be changed. Most importantly, the numerical steepest decent paths shall all be changed, and the stationary phase points and boundary resonance points are all different. As a result, the spacial relationships between stationary phase points and end points and numerical steepest decent paths and stokes line are all correspondingly changed. Therefore, the integral paths after adopting contour deformation technique are also changed. The NSDP method for the concave surface is proposed to compute the canonical PO type integral in a frequency-independent workload.

##### 3.2. Reduction of the Double Integral into Line Integral

The amplitude term of (10), , has the expression

Then, after integration by part, we obtain

Thenwith

In (13), there are two terms defined on the edge and another two terms defined on the edge , and the term defined on is easy to calculate. Therefore, is divided into and . They can be expressed as

Hence, we can rewrite as

In (15), both and are highly oscillatory integrals, but has the form of complementary error function and can be easily calculated with the aid of complementary error function. , however, can not be calculated analytically and the computational cost increases dramatically with the larger *k*. In the following, numerical steepest descent path method is illustrated to calculate .

##### 3.3. Calculating the Highly Oscillatory Integral

The contour deformation path technique is adopted to calculate . After substituting , we could obtain the phase function term of the integrand for as

The proposed steepest descent path should satisfy the following three conditions, :(a), which means that the NSDPs start from the end points , (b), where is a constant(c)

In this sense, the exponential term of the phase function can be derived as

It can be seen that term is a certain imaginary number and is a real term that decreases exponentially when goes large. Therefore, the exponential of phase function is clearly no longer oscillatory. The contrast between Figures 2(c)–2(d) and Figures 2(e)–2(f) illustrates the effect of our NSDP method.

From the definition of NSDPs, we could obtain steepest descent path for end points as , :wherewhere is the stationary phase point (SPP), which satisfies . The phase behavior of the SPP differs from the end points. In this sense, we define a contour path function for SPP as :

We could also obtain NSDP for as follows:

By adopting the NSDP method, can be written as

At this point, the behavior of complementary error function should be considered. has the following feature:where and are two slowly varying functions.

The above different behavior of on two different regions is called Stokes phenomenon. The line separating these two regions is called Stokes line, which can be expressed as

After substituting into (25), the Stokes line takes the following expression:

Considering the features of erfc in (23) can be further written aswhere

For integral , the original function of its integrand isand is from the constant “2” in (24), and its integral path is determined by the spacial relationship between Stokes line and NSDPs. Specifically, different cases of integral paths arise from relative locations of end points, SPP, and intersection point of Stokes line and *x* axis. Then can be calculated.

For , its specific expression comes from the relative values of end points , and SPP , which can be obtained by the following equations:

After substituting with the corresponding path function, we could obtain the expression for the terms in (30):

Here,

At this point, and are obtained, and can be calculated efficiently.

#### 4. An Example of Calculating PO Surface Integral by Adopting the NSDP Method

We consider a quadratic concave triangular patch , as shown in Figure 1(a). A flat triangular patch with three vertex points , , could be derived by projecting onto the - plane, as shown in Figure 1(b). And the three edges of are governed by ; the phase functions on these edges are , . We define , , and as edge 1, edge 2, and edge 3, respectively. As can be calculated analytically, in the following, we will discuss on three different edges of the triangular patch, respectively. First, taking the integral on edge 1 as an example, its steepest descent paths and Stokes line are shown in Figure 3(a).

**(a)**

**(b)**

**(c)**

For edge 1, which is , , , and correspond to , , and in the last section, respectively. Then, in (16) takes the following expression:

As in Figure 3(a), ; according to (30), can be written aswhere and have the similar forms as and given in (31).

As for , is the intersection point of Stokes line and steepest descent paths. Considering the integral path is in the domain satisfying , which is the domain above Stokes line, we can obtain as follows:where and have the similar forms as given in (30). At this point, for edge 1 is obtained.

Similarly, for edge 2, (, , and ), we have

To be noticed, as shown in Figure 3(b), . Therefore, according to (29), contains three terms, which are , , and , respectively.

And for edge 3, (, , and ), we have

Here,

Based on (42)–(44), we could derive the closed form formulas for defined on quadratic concave triangular patches:where is canceled during the NSDP technique process. We can also separate the contributions of three critical points, which are the stationary phase point (SPP), the resonance points (RP), and the vertex points (VP):

Hence, the PO surface integral can be expressed by the high frequency critical point contributions as follows:

#### 5. Numerical Results

##### 5.1. Calculating the Assembled Concave Quadratic Triangular Patches

We consider a concave quadratic patch (see Figure 4):with four vertex points: , , , and . The incident wave direction with the amplitude are considered when we treat a bistatic situation. The far field observation point lies along the direction . The quadrilateral domain in Figure 4(b) is decomposed into two triangular domains, and , respectively.

**(a)**

**(b)**

On invoking the affine transformation, we obtain the transformed new quadratic patch. The new vertex points for the transformed new patch are .

The concave quadratic patch can be calculated through two assembled triangular patches based on the NSDP method introduced above. The computational results and the corresponding consumed CPU time with increasing wavenumber are also calculated by the high frequency asymptotic (HFA) method and the brute force (BF). Specifically, from Figure 5(a), it can be inferred that the computational results of the NSDP method agree well with those from BF method and the HFA method. Furthermore, Figure 5(b) demonstrates that the NSDP method is able to calculate the scattered fields from concave surfaces in a frequency-independent manner.

**(a)**

**(b)**

##### 5.2. The Efficiency of the NSDP Method

We consider another concave PEC patch:to benchmark the efficiencies of our proposed NSDP method on calculating the PO scattered fields.

The projection of the surface on the plane is a quadrilateral domain, as shown in Figure 6, with four vertex points: , , , and . After the affine transformation, the quadrilateral domain in Figure 6(a) is transformed into another quadrilateral domain in Figure 6(b). We set the parameters as follows: the wave frequency and the incident wave propagates along direction. Besides the NSDP method, both the BF method and the HFA method are applied to calculate the PO scattered fields to illustrate the efficiency of the NSDP method.

**(a)**

**(b)**

Through the BF method verification, we could compare the results generated by the NSDP method and the HFA method. The computational results of the PO scattered field and the relative error generated by the NSDP and the HFA methods are shown in Table 1 and Table 2, respectively. Meanwhile, the CPU time consumed by these methods is given in Table 3. Figure 7(a) illustrates the RCS values of the PO scattered field produced by both methods, which match well with results from the BF method. Moreover, Figure 7(b) indicates that the accuracy can be significantly improved by around two digits through the NSDP method. Figure 7(c) demonstrates the CPU time consumed by both the NSDP method and the HFA method. Clearly, the computational workload for the NSDP method is relatively low and frequency-independent.

**(a)**

**(b)**

**(c)**

In summary, the proposed NSDP method can calculate the PO scattered fields from the quadratic concave surface in the frequency-independent computational workload and error-controllable accuracy.

#### 6. Conclusion

In this work, the NSDP method is proposed to compute the highly oscillatory PO scattered fields from concave quadratic surfaces, including both the monostatic and the bistatic cases. The considered concave surface is discretized into several triangular patches. The Lagrange interpolation polynomial approximation and the affine transformation are used to approximate the amplitude function and phase function of the PO type integral. By adopting the NSDP method to obtain corresponding NSDPs, highly oscillatory PO integrand defined on the assembled triangular patches is converted to exponentially decay integrand. Numerical examples indicate the comparison of the RCS values, the CPU time, and the relative error by the NSDP method, the HFA method, and the BF method. We make the conclusion that the proposed NSDP method for calculating the PO scattered fields is error-controllable and frequency-independent for the concave quadratic patches.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported in part by National Natural Science Foundation of China, No. 11671099; in part by National Science Fund for Distinguished Young Scholars, China, No. 61425010; in part by Science Challenge Project under Grant TZ2018002; in part by Open Research Fund of State Key Laboratory of Pulsed Power Laser Technology, China, No. SKL2018KF08; and in part by the Innovation Development Fund of CAEP, China, No. CX2019034.