Although FMM with ABC can improve the convergence rate, the parameter *h* is still a core factor affecting the value of convergence, as will be demonstrated in the next section. Moreover, the criterion for determining the appropriate *h* remains to be studied in detail. So far, we have implemented a computer program based on FMM with ABC for calculating the scattering characteristics of a graphene grating as well as the electric-field ((E_x)) distribution on the graphene grating surface. Let us skip ahead to the numerical result concerning a graphene-strip grating normally incident by a TM-polarized plane wave. Figure 5 shows the distribution of (E_x) (red dotted curve obtained using FMM with ABC) over the graphene-strip ((xin [0,20])) and slit ((xin [20,70])) regions. It is obvious to see jump discontinuities occurring around (x=0) and (x=20mu m); wherein (E_x) (or (J_x/sigma _g)) vanishes close to strip edges. In fact, induced current vanishes at strip edges can be explained physically as follows.

Referring to Fig. 1, at (x=w_g^+) in the slit region, magnetic field component (H_y) must be continuous across the interface due to (J_x(x=w_g^+,z=0)=0), namely, (H_y(x=w_g^+,z=0^-)=H_y(x=w_g^+,z=0^+)). Moreover, because of (H_y(x=w_g^+,z=0^-)=H_y(x=w_g^-,z=0^-)) in the lower medium and (H_y(x=w_g^+,z=0^+)=H_y(x=w_g^-,z=0^+)) in the upper medium, we have (H_y(x=w_g^-,z=0^-)=H_y(x=w_g^-,z=0^+)), resulting in the vanished induced current density (J_x(x=w_g^-,z=0)) obtained via Eq. (9); therefore we have (E_x(x=w_g^-,z=0)=0), that is, vanishing (E_x) at the edge. Additionally, the vanishing current at metal-strip edges was reported in the research of a metal-strip grating illuminated by a plane wave in TM polarization^{19}, but unfortunately the distribution of (E_x) in the slit region was not shown in that paper. Regarding the electric field in the slit region, an exponential growth of (E_x) around the slit edges can be observed (red dotted curve) in Fig. 5. Alternatively, as is well know, TM-polarized electric field near the edge of a thin sheet is proportional to (rho ^{-1/2}) and becomes singular as (rho ) approaches zero^{21}, where (rho ) is defined as the radius (in the polar coordinate system) with its original point locating at the edge.

Nevertheless, the Gibbs phenomenon taking place near the discontinuities reveals that the customary global basis of Floquet-Fourier series is inappropriate for expanding the field directly on the graphene grating surface. In view of that, it is essential to construct the local basis functions inherently satisfying the field nature in respective regions. Furthermore, the criterion for choosing basis functions contains: (1) to use only a few basis functions to approach the correct solution, and (2) to have closed forms in the overlap integral between the local basis functions and the space harmonic in Eq. (5).

More specifically, the (sin )-based local basis function vanishing at its both ends for any harmonic order is used to expand (E_x) over the graphene strip, which is given as

$$begin{aligned} g_n(x)=sqrt{frac{2}{w_g}}sin frac{npi (x-x^{(g)}_1)}{w_g}, end{aligned}$$

(23)

where (w_g=x^{(g)}_2-x^{(g)}_1); the graphene strip belongs to the region of ([x^{(g)}_1, x^{(g)}_2]); index *n* is ranging from 1 to (N_g).

On the other hand, in the slit region, we have the singular basis functions with singularities at its two edges, which are commonly used to approximate the current parallel to the edges in a micro-strip line^{22}. They are expressed as follows:

$$begin{aligned} s_n(x)=sqrt{frac{gamma _n}{w_s}}frac{cos frac{npi (x-x^{(s)}_1)}{w_s}}{sqrt{(w_s/2)^2-(x-x^{(s)}_c)^2}}, end{aligned}$$

(24)

where (x^{(s)}_c=(x^{(s)}_1+x^{(s)}_2)/2) and (w_s=x^{(s)}_2-x^{(s)}_1); the slit is in the region of ([x^{(s)}_1, x^{(s)}_2]). Parameter (gamma _n=2) for (n=0) and (gamma _n=1) for (nne 0); index *n* runs from 0 to (N_s-1). Notably, the denominator in (s_n(x)) approximates (sqrt{w_s}sqrt{rho }) at the strip edge where (x= w_g+rho ) ((x_1^{(s)}=w_g) and (x_2^{(s)}=d_x)), which confirms the electric-field edge condition described previously. Notably, Eqs. (23) and (24) both are expressed in the general form for easy extension to the case with multiple graphene strips in a period.

In a unit cell on the graphene grating surface (at (z=0)), (E_x(x)) can be written as

$$begin{aligned} E_x(x)= left{ begin{array}{rcl} sum limits _{n=1}^{N_g} p_n g_n(x) &{} text{ for } &{} xin text{ graphene } \ sum limits _{n=0}^{N_s-1} q_n s_n(x) &{} text{ for } &{} xin text{ slit } end{array}. right. end{aligned}$$

(25)

Parameters (N_g) and (N_s) represent the number of basis in the graphene strip and slit regions, respectively.

Due to electromagnetic boundary condition, (E_x) must be continuous at the interface between graphene grating and uniform medium at (z=0). Therefore, equality of Eqs. (4) and (25) gives

$$begin{aligned} sum _{n=-N}^{n=+N}V_n(0)psi _n(x)=left{ begin{array}{rcl} sum limits _{n=1}^{N_g} p_n g_n(x) &{} text{ for } &{} xin text{ graphene } \ sum limits _{n=0}^{N_s-1} q_n s_n(x) &{} text{ for } &{} xin text{ slit } end{array}right. . end{aligned}$$

(26)

By multiplying the complex conjugate of (psi _m(x)) on both sides of Eq. (26) and taking integration over one period, we obtain

$$begin{aligned} sum _{n=-N}^{n=+N}V_n(0)langle psi _m(x)^dagger | psi _n(x)rangle =sum _{n=1}^{N_g}p_n langle psi _m(x)^dagger | g_n(x)rangle +sum _{n=0}^{N_s-1} q_n langle psi _m(x)^dagger | s_n(x)rangle , end{aligned}$$

(27)

where integer *m* is ranging from –*N* to +*N*. The notation of (langle a(x) | b(x) rangle =int _{0}^{d_x}a(x)b(x) dx) is defined as the overlap integral of functions *a*(*x*) and *b*(*x*) in the range of ([0,d_x]). By the orthogonality of (psi _n(x)) in Eq. (6) and the closed form solutions of overlap integral, the above equation becomes

$$begin{aligned} V_m(0)=sum _{n=1}^{N_g}p_n G_{mn}+sum _{n=0}^{N_s-1}q_n S_{mn}, end{aligned}$$

(28)

with

$$begin{aligned} G_{mn}= & {} langle psi _m(x)^dagger | g_n(x)rangle = frac{-j}{sqrt{2}}sqrt{frac{w_g}{d_x}} e^{jk_{xm} x^{(g)}_c} cdot left[ e^{jnpi /2} sinc(alpha _{mn}^{+} w_g/2)+ e^{-jnpi /2} sinc(alpha _{mn}^{-} w_g/2) right] , end{aligned}$$

(29)

$$begin{aligned} S_{mn}= & {} langle psi _m(x)^dagger | s_n(x)rangle = frac{pi }{2}sqrt{frac{gamma _n}{w_s d_x}}e^{jk_{xm}x^{(s)}_c} cdot left[ e^{jnpi /2} J_o(beta _{mn}^{+}w_s/2)+ e^{-jnpi /2} J_o(beta _{mn}^{-} w_s/2) right] , end{aligned}$$

(30)

where (beta _{mn}^{pm }=k_{xm}pm npi /w_s) and (alpha _{mn}^{pm }=k_{xm} pm npi /w_g); function (J_o(cdot )) is the zero order Bessel function of the first kind; function *sinc*(*x*) is the unnormalized sinc function defined as (sinc(x)=sin(x)/x). Parameter *n* is the index of the local basis function.

Equation 28 can be rewritten as a matrix-vector form:

$$begin{aligned} {underline{V}}(0) = [[G]]{underline{p}}+[[S]]{underline{q}}= begin{bmatrix} {[}[G]]&[[S]] end{bmatrix} begin{bmatrix} {underline{p}} \ {underline{q}} end{bmatrix} . end{aligned}$$

(31)

Vector ({underline{p}}) is a (N_g)-by-1 column vector with its (n{text{th}}) element (p_n); (q_n) is the (n{text{th}}) element in column vector ({underline{q}}) of size (N_s)-by-1. Here, the sub-matrix [[*G*]] and [[*S*]] have the size (N_{tot})-by-(N_g) and (N_{tot})-by-(N_s), respectively. Specifically, (N_{tot}), (N_g) and (N_s) satisfy the relationship: (N_{tot}=N_g+N_s). Moreover, the ratio between (N_g) and (N_s) equals to the ratio of (w_g) to (w_s)^{5, 19}, therefore, we have (N_g=round[N_{tot}cdot w_g/(w_g+w_s)]) and (N_s=N_{tot}-N_g); the operator *round*[.] rounds a real number towards the nearest integer. In doing so, we have a square matrix (left[ [G]] [[S] right] ) of size (N_{tot})-by-(N_{tot}).

Substituting of Eqs. (3) and (25) into Eq. (9), we obtain

$$begin{aligned} sum _{n=-N}^{n=+N}[I_n(0^-)-I_n(0^+)]psi _n(x) =sum _{n=1}^{N_g} sigma _g p_n g_n(x). end{aligned}$$

(32)

Multiplying (psi _m^dagger (x)) on both sides of Eq. (32) and taking the integration over one period along the *x*-axis, one obtains

$$begin{aligned} sum _{n=-N}^{n=+N}[I_n(0^-)-I_n(0^+)]langle psi _m(x)^dagger | psi _n(x)rangle =sum _{n=1}^{N_g} sigma _g p_n langle psi _m(x)^dagger | g_n(x)rangle . end{aligned}$$

(33)

By invoking orthogonality and Eq. (29), the system of linear equations in Eq. (33) can be expressed in terms of matrix-vector form. One obtains

$$begin{aligned} {underline{I}}(0^{-})-{underline{I}}(0^{+})= begin{bmatrix} sigma _g[[G]]&[[0]] end{bmatrix} begin{bmatrix} {underline{p}} \ {underline{q}} end{bmatrix} , end{aligned}$$

(34)

where [[0]] is a null matrix of size (N_{tot})-by-(N_s), and matrix (left[ sigma _g[[G]] [[0]]right] ) is a square matrix of size (N_{tot})-by-(N_{tot}).

After performing matrix operations with Eqs. (31) and (34), we obtain the Eq. (12) with a new admittance matrix of size (N_{tot})-by-(N_{tot}) given below

$$begin{aligned} {[}[Y_g]]= begin{bmatrix} sigma _g[[G]]&[[0]] end{bmatrix} begin{bmatrix} {[}[G]]&[[S]] end{bmatrix} ^{-1}, end{aligned}$$

(35)

where the new admittance matrix in Eq. (35) is a square matrix of size (N_{tot})-by-(N_{tot}).

Here, we obtain a totally different admittance matrix while the input-output relation in Eq. (12) remains the same. The same procedure in FMM can be applied to calculate the reflect and transmit amplitudes of each diffraction order. Once the voltage ({underline{V}}(0)=[[T]]{underline{a}}) is obtained, the expansion coefficients of local basis functions, ({underline{p}}) and ({underline{q}}) in Eq. (31) can be readily determined, as well as the distribution of (E_x) on the graphene strip and slit.

Furthermore, the power dissipated on the graphene strips array can be directly determined by

$$begin{aligned} P_{abs.}=frac{1}{2}{text {Re}}left[ sum _{n=1}^{n=N_g} sigma _g mid p_n mid ^2 right] . end{aligned}$$

(36)