The Orthogonal Expansion in Time-Domain Method for Solving Maxwell Equations Using Paralleling-in-Order Scheme

Abstract

The orthogonal expansion in time-domain method is a new kind of unconditionally stable finite-difference time-domain (FDTD) method for solving the Maxwell equation efficiently. Generally, it can be implemented by two schemes: marching-on-in-order and paralleling-in-order, which, respectively, use weighted Laguerre polynomials and associated Hermite functions as temporal expansions and testing functions. This chapter summarized paralleling-in-order-based FDTD method using associated Hermite functions and Legendre polynomials. And a comparison from theoretical analysis to numerical examples is shown. The LD integral transfer matrix can be considered as a “dual” transformation for AH differential matrix, which gives a possible way to find more potential orthogonal basis function to implement a paralleling-in-order scheme. In addition, the differences with these two orthogonal functions are also analyzed. From the numerical results, we can see their agreements in some general cases while differing in some cases such as shielding analysis with the long-time response requirement.

1. Introduction

To overcome the numerical stability constraints of conventional finite-difference time-domain (FDTD) method [1, 2], many unconditionally stable methods to reduce or eliminate requirements of the stability condition have been proposed and developed, such as alternating-direction implicit method [2, 3] and locally one-dimensional schemes [3], explicit and unconditionally stable FDTD method [4], and orthogonal expansions in time domain [5, 6, 7, 8]. For the orthogonal expansions schemes, field-versus-time variations in the FDTD space lattice are expanded using an appropriate set of orthogonal temporal basis and testing functions, such as weighted Laguerre polynomials (WLP) and associated Hermite (AH) functions, which leads to two different solution schemes: marching-on-in-order and paralleling-in-order, respectively. Both of them appear to be promising according to the reported work where the computational time can be reduced to at least 10% of the conventional FDTD scheme [1]. Recently, the Legendre (LD) polynomials are explored as another possible orthogonal expansion incorporated with FDTD to form a paralleling-in-order-based unconditionally stable FDTD method. Based on it, in this chapter, we made a comparison investigation for these two new methods, which are AH FDTD method and LD FDTD method, especially focused on their differences. Through a numerical example, we validate their effectiveness when compared with the conventional FDTD method and summarized the characteristics of the two methods.

2. Formulation for paralleling-in-order scheme: AH and LD functions

2.1 2D Maxwell’s equations in time domain

The 2D time-domain Maxwell’s equations with the TEz wave case in lossy medium are considered:

$$\epsilon \frac{\partial E_x\left(r,t\right)}{\partial t}+{\sigma }_e{E}_x\left(r,t\right)=\frac{\partial H_z\left(r,t\right)}{\partial y}-J_x\left(r,t\right)$$ E1

$$\mu \frac{\partial H_z\left(r,t\right)}{\partial t}+{\sigma }_m{H}_z\left(r,t\right)=\frac{\partial E_x\left(r,t\right)}{\partial y}-\frac{\partial E_y\left(r,t\right)}{\partial x}-M_z\left(r,t\right)$$ E2

$$\epsilon \frac{\partial E_y\left(r,t\right)}{\partial t}+{\sigma }_e{E}_y\left(r,t\right)=-\frac{\partial H_z\left(r,t\right)}{\partial x}-J_y\left(r,t\right)$$ E3

where $\epsilon$ , $\mu$ , ${\sigma }_e$ , and ${\sigma }_m$ are the permittivity, the permeability, the electric conductivity, and the magnetic loss of the medium, respectively. $E_{\xi }\left(r,t\right)$ and $J_{\xi }\left(r,t\right)$ ( $\xi =x,y$ ) are the electric field component and the electric current densities, respectively. $H_z\left(r,t\right)$ and $M_z\left(r,t\right)$ are the magnetic field component and magnetic current densities, respectively.

2.2 The differential and integral transfer matrices to deal with the partial differential term in Maxwell’s equations

2.2.1 The associated Hermite function

Associated Hermite function is defined as

$$\left\{{\varphi }_n\left(t\right)={\left(2^{n}n!{\pi }^{1/2}\right)}^{-1/2}{e}^{-t^2/2}{H}_n\left(t\right)\right\},\left(\mathrm{n}=0,1\dots \right)$$ E4

where $H_n\left(t\right)={\left(-1\right)}^n{e}^{t^2}\frac{d^n}{{\mathit{dt}}^n}\left(e^{-t^2}\right)$ is Hermite polynomials. Although it is not causal, it can be transformed into causal form by virtue of a proper translating and scaling parameters and then used to span the causal electromagnetic responses. The transformed basis function is $\left\{{\overline{\varphi }}_n\left(\tilde{t}\right)={\left(2^{n}n!{\mathit{\sigma \pi }}^{1/2}\right)}^{-1/2}{e}^{-\tilde{t}{}^2/2}{H}_n\left(\tilde{t}\right)\right\}$ , where transformed time variable $\tilde{t}=\left(t-T_f\right)/\sigma$ . And $T_f$ is a translating parameter and $\sigma$ is a scaling parameter. By controlling these two parameters, the time-frequency support of the AH functions $\left\{{\overline{\varphi }}_n\left(\tilde{t}\right)\right\}$ space can be changed flexibly. So, arbitrary locally time-supported functions can be spanned by these transformed basis functions, including the causal electromagnetic responses.

From [7], if a causal function $u\left(r,t\right)$ , such as the electric or magnetic field function, can be expanded by

$$u\left(r,t\right)=\sum _{n=0}^{\infty }{u}_n\left(r\right){\overline{\varphi }}_n\left(\tilde{t}\right)$$ E5

we can deduce the first derivative of $u\left(x,t\right)$ with respect to

$$\frac{\partial }{\partial t}u\left(r,t\right)=\frac{1}{\sigma }\sum _{n=0}^{\infty }\left(u_{n+1}\left(r\right)\sqrt{\frac{n+1}{2}}-u_{n-1}\left(r\right)\sqrt{\frac{n}{2}}\right){\overline{\varphi }}_n\left(\tilde{t}\right)$$ E6

Then, the Q-tuple AH domain coefficients for $u\left(r,t\right)$ and $\dot{u}\left(r,t\right)$ from (5) and (6) can be obtained as $U={\left[U^0\cdots U^{Q-1}\right]}^T$ and $\dot{U}={\left[{\dot{U}}^0\cdots {\dot{U}}^{Q-1}\right]}^T$ . And, we can readily obtain the relationship between $U$ and $\dot{U}$ as

$$\dot{U}=\mathit{\alpha U}$$ E7

where

$$\alpha =\frac{\sqrt{2}}{2\lambda }{\left[\begin{array}{l}\sqrt{1}\\ -\sqrt{1}\sqrt{2}\\ -\sqrt{2}\ddots \\ \ddots \sqrt{Q-1}\\ -\sqrt{Q-1}\end{array}\right]}_{Q\times Q}$$ E8

By using (8), the partial differential term in Maxwell?s equations can readily be dealt with, and finally, a five-point banded matrix equation for Hz component can be obtained [9].

2.2.2 The associated Legendre polynomial

We expand all the temporal quantities in terms of the associated Legendre polynomial given by [10]:

$$P_q\left(t\right)=\sqrt{\frac{2q+1}{2}}{L}_q\left(2\frac{t}{l}-1\right),t\in \left[0,l\right]$$ E9

where l is the time support for analyzing a causal response and Lq is the Legendre polynomial with order q, which are orthogonal in the interval [−1,1] satisfying the following recurrence relation:

$$L_{q+1}\left(t\right)=\frac{2q+1}{q+1}{\mathit{tL}}_q\left(t\right)-\frac{q}{q+1}{L}_{q-1}\left(t\right),$$ E10

and $L_0\left(t\right)=0$ , $L_1\left(t\right)=t$ . Given a time-support field function $u\left(r,t\right)$ , it can be expanded by (9) as

$$u\left(r,t\right)=\sum _{q=0}^{\infty }{u}_q\left(r\right)P_q\left(t\right)$$ E11

where $u_q\left(r\right)$ is the q-th expanding coefficients, and it can be calculated by

$$u_q\left(r\right)=\underset{-\infty }{\overset{+\infty }{\int }}u\left(r,t\right)P_q\left(t\right)\mathit{dt}$$ E12

From the intrinsic features of Legendre function, the differential relationship can be described as

$$P_q\left(t\right)=\frac{1}{2\sqrt{\left(2q+3\right)\left(2q+1\right)}}{P}_{q+1}^{\text{'}}\left(t\right)-\frac{1}{2\sqrt{\left(2q+1\right)\left(2q-1\right)}}{P}_{q-1}^{\text{'}}\left(t\right)$$ E13

If the field derivative of $u\left(r,t\right)$ to t is expanded as

$$u^{\text{'}}\left(r,t\right)=\sum _{q=0}^{\infty }{u}_q^{\left(1\right)}\left(r\right)P_q\left(t\right)$$ E14

where $u_q^{\left(1\right)}\left(r\right)$ is q-th expanding coefficients for $u^{\text{'}}\left(r,t\right)$ , then incorporated with (13), it can be deduced as

$$u^{\text{'}}\left(r,t\right)={\left(\sum _{q=0}^{\infty }\left(\frac{l}{2\sqrt{\left(2q+1\right)\left(2q-1\right)}}{u}_{q-1}^{\left(1\right)}\left(r\right)-\frac{l}{2\sqrt{\left(2q+3\right)\left(2q+1\right)}}{u}_{q+1}^{\left(1\right)}\left(r\right)\right)P_q\left(t\right)\right)}^{\text{'}}$$ E15

Connecting (15) and (11), we can get

$$u_q\left(r\right)=\frac{l}{2\sqrt{\left(2q+1\right)\left(2q-1\right)}}{u}_{q-1}^{\left(1\right)}\left(r\right)-\frac{l}{2\sqrt{\left(2q+3\right)\left(2q+1\right)}}{u}_{q+1}^{\left(1\right)}\left(r\right)$$ E16

When assembling ${\left\{u_q\left(r\right)\right\}}_{q=0,1\cdots Q-1}$ as a Q-tuple $U$ and ${\left\{u_q^{\left(1\right)}\left(r\right)\right\}}_{q=0,1\cdots Q-1}$ as $U^{\left(1\right)}$ , a matrix-multiply relationship can be obtained from (16) as the following:

$$U={\alpha }_L{U}^{\left(1\right)}$$ E17

where ${\alpha }_L$ is integral matrix.

$${\alpha }_L=\frac{l}{2}{\left[\begin{array}{l}‐1/\sqrt{1\cdot 3}\\ 1/\sqrt{1\cdot 3}‐1/\sqrt{3\cdot 5}\\ 1/\sqrt{3\cdot 5}\ddots \\ \ddots ‐1/\sqrt{\left(2Q-3\right)\left(2Q-1\right)}\\ 1/\sqrt{\left(2Q-3\right)\left(2Q-1\right)}\end{array}\right]}_{Q\times Q}$$ E18

Alternatively, Eq. (17) can be rewritten as.

$$U^{\left(1\right)}={\alpha }_L^{-1}U$$ E19

2.3 From time domain to orthogonal domain and reconstruction

When the differential or integral transfer matrices are obtained, the time-domain Maxwell equation can be transformed directly into AH or LD domain. Here, let us set LD as an example to illustrate the later formulation.

Similar to the paralleling-in-order-based AH FDTD method, we can apply a Q-tuple-domain transformation for LD FDTD method to (1)–(3) and discretize them as the following:

$${\alpha }_{e\left(i,j\right)}{\left.E_x\right|}_{i,j}=\left({\left.H_z\right|}_{i,j}-{\left.H_z\right|}_{i,j-1}\right)/\mathrm{\Delta }{\overline{y}}_j-{\left.J_x\right|}_{i,j}$$ E20

$${\alpha }_{e\left(i,j\right)}{\left.E_y\right|}_{i,j}=-\left({\left.H_z\right|}_{i,j}-{\left.H_z\right|}_{i-1,j}\right)/\mathrm{\Delta }{\overline{x}}_i-{\left.J_y\right|}_{i,j}$$ E21

$${\left.{\alpha }_{m\left(i,j\right)}{H}_z\right|}_{i,j}=\left({\left.E_x\right|}_{i,j+1}-{\left.E_x\right|}_{i,j}\right)/\mathrm{\Delta }{y}_j-\left({\left.E_y\right|}_{i+1,j}-{\left.E_y\right|}_{i,j}\right)/\mathrm{\Delta }{x}_i-{\left.M_z\right|}_{i,j}$$ E22

where

$${\alpha }_{e\left(i,j\right)}={\left.\epsilon \right|}_{i,j}{\alpha }_L^{-1}+{\left.{\sigma }_e\right|}_{i,j}I$$ E23

$${\alpha }_{m\left(i,j\right)}={\left.{\mu }_m\right|}_{i,j}{\alpha }_L^{-1}+{\left.{\sigma }_m\right|}_{i,j}I$$ E24

where ${\left.E_x\right|}_{i,j}$ , ${\left.E_y\right|}_{i,j}$ , ${\left.H_z\right|}_{i,j}$ , ${\left.J_x\right|}_{i,j}$ , ${\left.J_y\right|}_{i,j}$ , and ${\left.M_z\right|}_{i,j}$ are Q-tuple representations of fields and sources, respectively. And, $I$ is the Q-dimensional identity matrix. By assembling (20)–(22) and eliminating the electric field components, a five-diagonal banded matrix equation for Hz component can be obtained:

$$a_{l\left(i,j\right)}{\left.H_z\right|}_{i-1,j}+a_{r\left(i+1,j\right)}{\left.H_z\right|}_{i+1,j}+a_{m\left(i,j\right)}{\left.H_z\right|}_{i,j}+a_{d\left(i,j\right)}{\left.H_z\right|}_{i,j-1}+a_{u\left(i,j+1\right)}{\left.H_z\right|}_{i,j+1}=b_{i,j}$$ E25

where

$$a_{u\left(i,j+1\right)}=-{\alpha }_{e\left(i,j+1\right)}^{-1}/\mathrm{\Delta }{\overline{y}}_{j+1}/\mathrm{\Delta }{y}_j$$ E26

$$a_{d\left(i,j\right)}=-{\alpha }_{e\left(i,j\right)}^{-1}/\mathrm{\Delta }{\overline{y}}_j/\mathrm{\Delta }{y}_j$$ E27

$$a_{l\left(i,j\right)}=-{\alpha }_{e\left(i,j\right)}^{-1}/\mathrm{\Delta }{\overline{x}}_i/\mathrm{\Delta }{x}_i$$ E28

$$a_{r\left(i+1,j\right)}=-{\alpha }_{e\left(i+1,j\right)}^{-1}/\mathrm{\Delta }{\overline{x}}_{i+1}/\mathrm{\Delta }{x}_i$$ E29

$$a_{m\left(i,j\right)}=-\left(a_{r\left(i+1,j\right)}+a_{l\left(i,j\right)}+a_{u\left(i,j+1\right)}+a_{d\left(i,j\right)}+{\alpha }_{m\left(i,j\right)}\right)$$ E30

$$b_{i,j}=-\left({\left.{\alpha }_{m\left(i,j+1\right)}^{-1}{J}_x\right|}_{i,j+1}-{\left.{\alpha }_{m\left(i,j\right)}^{-1}{J}_x\right|}_{i,j}\right)/\mathrm{\Delta }y+\left({\left.{\alpha }_{m\left(i+1,j\right)}^{-1}{J}_y\right|}_{i+1,j}-{\left.{\alpha }_{m\left(i,j\right)}^{-1}{J}_y\right|}_{i,j}\right)/\mathrm{\Delta }x{\left.-M_z\right|}_{i,j}$$ E31

By using eigenvalue transformation from ${\alpha }_{L}X=\mathit{XV}$ , where $X$ and $V$ are the eigenvector matrix and diagonal matrix composed of eigenvalues $\left\{{\lambda }_q\right\}$ , respectively, Eq. (25) can be changed to the paralleling-in-order solution. For the q-th decoupled equation, we have

$$A\left(1/{\lambda }_q\right){\left.H_z^{\ast }\right|}^q={\left.b^{\ast }\right|}^q$$ E32

where $A\left(·\right)$ is a banded sparse matrix, with the similar form as from AH FDTD method, and ${\left.b^{\ast }\right|}^q$ is the transformed variables from ${\left.b\right|}_{i,j}=X{\left.b^{\ast }\right|}_{i,j}$ . Finally, we can obtain a paralleling-in-order scheme to calculate all of the expanding coefficients of electromagnetic fields, and then the time-domain responses can be reconstructed from (11).

3. Comparison for the two methods

The above formula can be regarded and classified as a uniform OF differential transfer matrix transformation. Therefore, as long as the LD differential matrix is replaced by the AH domain differential transfer matrix, the FDTD algorithm based on the LD orthogonal basis function, LD FDTD, including the parallel solution AH FDTD algorithm [9], and the alternate direction efficient calculation [11] can be easily realized. The implementation of the program only requires a simple modification.

Table 1 gives a comparison of the relevant properties of the LD FDTD method and the AH FDTD method. It can be seen that the two methods can be considered as a “dual” system, because the AH differential matrix is the basic element of the AH FDTD method and the LD integration matrix is also the basic element of the LD FDTD method. This gives us a revelation that is it possible that any orthogonal basis function can construct a differential or integral transfer matrix and then easily implement a paralleling-in-order scheme similar like AH FDTD algorithm? The answer might be NOT. Such as the Laguerre FDTD method, as introduced before, cannot be calculated in parallel. However, it is undeniable that there may be more basis functions that can implement the paralleling-in-order scheme. If any, we can collectively call these methods as the AH series unconditionally stable FDTD method.

AH FDTD	LD FDTD
Differential transfer matrix	Integral transfer matrix
$U^{\left(1\right)}=\mathit{\alpha U}$	$U={\alpha }_L{U}^{\left(1\right)}$
${\alpha }_{\left(l\right)}=\frac{\sqrt{2}}{2l}{\left[\begin{array}{l}\sqrt{1}\\ -\sqrt{1}\sqrt{2}\\ -\sqrt{2}\ddots \\ \ddots \sqrt{Q-1}\\ -\sqrt{Q-1}\end{array}\right]}_{Q\times Q}$	${\alpha }_{L\left(l\right)}=\frac{l}{2}{\left[\begin{array}{l}‐1/\sqrt{1\cdot 3}\\ 1/\sqrt{1\cdot 3}‐1/\sqrt{3\cdot 5}\\ 1/\sqrt{3\cdot 5}\ddots \\ \ddots ‐1/\sqrt{\left(2Q-3\right)\left(2Q-1\right)}\\ 1/\sqrt{\left(2Q-3\right)\left(2Q-1\right)}\end{array}\right]}_{Q\times Q}$
$\frac{\partial }{\partial t}\leftarrow \alpha \to \mathit{j\omega }$	$\int \mathit{dt}\leftarrow {\alpha }_L\to \frac{1}{\mathit{j\omega }}$
$\left\{\begin{array}{l}{T}_Q\approx 2l\left(\sqrt{\mathit{\pi Q}/1.7}+1.8\right)\\ F_Q\approx \frac{\sqrt{\mathit{\pi Q}/1.7}+1.8}{2\mathit{\pi l}}\end{array}\right.\to \left(l,Q\right)$	Scale factor l = TQ; Finite order of Q
With time-frequency Homomorphism	Without time-frequency; Homomorphism
Antisymmetry; Eigenvalue conjugate symmetry	Antisymmetry; Eigenvalue conjugate symmetry
$A\left({\lambda }_q\right)H^q=J^q$	$A\left(1/{\lambda }_q\right)H^q=J^q$

Table 1.

LD comparison of LD FDTD method and AH FDTD method.

4. Numerical verification

4.1 An infinitely large lossy dielectric plate

As AH or LD FDTD method shares with almost the same program, a 1-D program is set for a general verification. Figure 1 shows the simulation results when a uniform plane wave penetrates an infinitely large lossy dielectric plate. The figure includes the electric field waveforms calculated by the AH FDTD method and the LD FDTD method and their relative errors with respect to the conventional FDTD method. It can be seen that the time-domain waveforms of both can be consistent with the results of the FDTD method and the relative errors are basically the same, only differing in the initial part. Therefore, in general, when the order of the two basic functions is the same and the parameters are selected reasonably, the accuracy is basically the same, and the efficiency is almost the same.

Figure 1.

Comparison of calculation results between AH FDTD method and HR FDTD method when simulating an infinitely large lossy dielectric plate. (a) Time-domain waveform. (b) Relative error.

4.2 An nonuniform parallel plate waveguide with a slot

However, the two methods also have the differences when simulating the long-time response applications, such as the example in [12]. The numerical example is set as a TEz wave propagation in a parallel plate waveguide, as shown in Figure 2. It is with a PEC slot of the thickness 0.2 mm and the distance 0.2 mm and a partly filled dielectric material of the thickness 0.8 mm with the dielectric medium parameters given as two cases: case I, ε = 11 ε₀, μ = μ₀, σ_e = 0.003 S/m, and σ_m = 0 Ω/m; case II, ε = 2 ε₀, μ = μ₀, σ_e = 30,000 S/m, and σ_m = 0 Ω/m. There are 140 × 8 uniform cells ( ${\mathrm{\Delta }}_x={\mathrm{\Delta }}_y$ =0.1 mm) in the computational domain. A Gaussian pulse sinusoidally modulated is used as the electric current source profile:

$$J_y\left(t\right)=\text{exp}\left(-{\left(\left(t-t_c\right)/t_d\right)}^2\right)\text{sin}\left(2\pi f_c\left(t-t_c\right)\right)$$ E33

Figure 2.

The geometry configuration for a 2D parallel plate waveguide with a PEC slot and a partly filled dielectric medium [12].

where $t_d$ = $1/\left(2f_c\right)$ , $t_c$ = 4 $t_d$ , and $f_c$ = 12 GHz. And the total simulation time is set as $l$ = 1.28 ns for case I and $l$ = 12.8 ns for case II; then it leads to the marching-in-on-time steps for N = 6000 and N = 60,000, respectively. And the number of orders for LD functions is chosen as 80 and 300, respectively, to obtain a good approximation of field components.

The Ey electric field responses at measurement point p1 and p2, located at the center of the slot and behind the medium, respectively, are calculated, which are both in agreement with the conventional FDTD method as shown in Figures 2 and 3. For comparison, the AH FDTD method is also used in these two cases. One can find the good results in Figure 3, but the errors come out in Figure 4 for AH FDTD method when the same number of orthogonal functions (Q = 80 for case I or 300 for case II) is used as LD FDTD method. However, when Q reaches 800, the results from AH FDTD method can achieve a comparable accuracy with the ones from LD FDTD method. One should note that for case II the waveform at point p2 has larger amplitude attenuation and longer delay than the result at point p1 due to the high dielectric medium located between them.

Figure 3.

The calculated results of transient electric field Ey for the case of I [12].

Figure 4.

The calculated results of transient electric field Ey for the case of II [12].

Tables 2 and 3 show the comparison of the computational resources. We can see that the simulation takes much more time for the FDTD method compared with proposed method, especially for the case of II, while the trade-off for the proposed method is that it consumes more memory than conventional FDTD method, which is similar to the AH FDTD method. In addition, from Table 3, we can find the advantages compared with AH FDTD method that the proposed method can use relative smaller memory storage and slightly fewer CPU times to get a readily results.

	Δt (ps)	Memory (MB)	CPU time (s)
FDTD (N = 6000)	0.21	1.8	2.97
AH FDTD (Q = 80)	21	2.9	1.32
LD FDTD (Q = 80)	21	2.9	1.32

Table 2.

The comparison of computational resources for the case of I [12].

	Δt (ps)	Memory (MB)	CPU time (s)
FDTD (N = 60,000)	0.21	1.8	30.8
AH FDTD (Q = 300)	21	11.8	1.55
AH FDTD (Q = 800)	21	28.9	1.95
LD FDTD (Q = 300)	21	11.8	1.55

Table 3.

The comparison of computational resources for the case of II [12].

5. Conclusions and future developments

The paralleling-in-order-based unconditionally stable FDTD methods are introduced using associated Hermite and Legendre polynomials in this chapter. The direct Q-tuple-domain transformation for time-domain Maxwell equation is guaranteed by using the integral matrix and differential matrix for Legendre function and associated Hermite functions that are introduced from the intrinsic integral or differential features for these orthogonal functions. Normally, the integral matrix of Legendre function can be considered as an inverse relationship from the differential operator, similar to the AH differential matrix. From this view, we can consider them as a uniform algorithm organized from the paralleling-in-order solution scheme. In addition, this chapter also detailed the different properties and the formula with these two methods theoretically and tested by numerical examples. Numerical examples for 1D and 2D cases validate their effectiveness and show LD FDTD with a better performance than AH FDTD method, in long-time simulation applications. In the next step, the more general paralleling-in-order scheme should be summarized, and then find or construct other possible orthogonal functions for their specific applications.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grants 61801217 and 51477183 and Natural Science Foundation of Jiangsu Province under Grant BK20180422. This support is gratefully acknowledged.