Applications of Tunable Mid-Infrared Plasmonic Square-Nanoring Resonator Based on Graphene Nanoribbon

In this work, different structures are designed based on graphene square-nanoring resonator (GSNR) and simulated by the three-dimensional finite-difference time-domain (3D-FDTD) method. Depending on the location and number of graphene nanoribbons (GNR), the proposed structures can be utilized as a band-pass filter, wavelength demultiplexer, or power splitter in the mid-infrared (MIR) wavelengths. The tunability of the suggested assemblies may be controlled simply by changing the dimensions and/or the chemical potential of the GSNRs. Benefiting from the nanoscale and ultra-compact GNRs, these structures can be proposed as basic blocks for optical computing and signal processing in the MIR wavelengths.


Introduction
Surface plasmon polaritons (SPPs), the electromagnetic waves propagating at the conductor-dielectric interfaces, are known as reliable and promising phenomena to overcome diffraction limit and realize devices in sub-wavelength and nanometer scales [1][2][3]. Diverse plasmonic devices have been designed and proposed based on controlling and manipulation of the propagation performance of SPPs at the metaldielectric surfaces. However, major barriers and difficulties, such as high loss and the lack of dynamic tunability, limit these plasmonic structures to be employed and implemented in practical applications [4][5][6][7]. Due to its unique electrical and optical properties, graphene has become an attractive material for researchers, in recent years, and can be used as a promising platform for photonic and plasmonic applications. This is due to its surface ability to support SPPs, known as graphene surface plasmon polaritons (GSPPs) [8][9][10], leading to significant progress observed in the field of this material plasmonics.
Graphene is a two-dimensional material that behaves like a single atom layer of metal since it has a negative relative permittivity in the mid-infrared (MIR) wavelength range [11,12]. The surface conductivity of graphene has a wide range of tunable characteristics, obtained only via the Fermi level adjustment of the material through electrostatic gating or chemical doping, which makes it widely used in tunable and high-performance optical devices [13]. The inclusion of graphene in the design and fabrication of integrated nanoplasmonic and nanophotonic devices brings many advantages compared to noble metals. To name a few [14][15][16]: (1) Ultra-thin thickness of graphene (much thinner than any metallic layer) makes it suitable to be used as waveguides in the supported wavelength ranges. (2) It has low ohmic losses at high-doped levels. (3) The large effective index of graphene structures in the MIR range, leading to high compression of the light and ultra-short local wavelengths within the material is a key to overcome the diffraction limit of the light compared to the traditional optical devices. In other words, SPP waves in graphene have wave vectors nearly 100 times larger than their free space counterparts, which makes the design of graphene-based optical devices easier, because the dimensions of the components are below the diffraction limit. (4) The conductivity of the material can be changed (tuning the chemical potential) by either doping or gate voltage (electrostatic gating). This makes the device to be set and calibrated without the need of re-fabricating the structure. (5) GSSPs have longer propagation lengths in the MIR region than metallic surface plasmon waves.
Over the two past decades, there has been extensive researches on a variety of graphene-based plasmonic devices such as modulators [17][18][19][20], switches [21,22], filters [23][24][25], power splitter [26,27], metamaterials [28][29][30][31], metasurfaces [32,33], as well as graphene-based plasmonic induced transparency (PIT) and reflection-induced (PIR) effects [34][35][36]. Most of these devices work based on graphene nanoribbons (GNRs), formed by reducing the width of the graphene flake. When the width of the graphene ribbon is reduced to less than a few tens of nanometers, the waveguide mode disappears and only the edge mode at the rims of the graphene strips can be excited vertically and laterally. In addition, the edge modes supported by GNRs have two distinct odd and even edge mode participation [37,38]. When the width of a GNR is a few tens of nanometers, the odd edge modes vanish, and we only have the even modes. This even edge mode has great transmission property, especially in the field of filter design. Another advantage of GNR in such devices is the relatively higher effective refractive index of the edge modes compared to waveguide modes, and consequently, a much smaller effective wavelength [39,40]. This makes the GNR-based plasmonic device ultra-compact with a small footprint.
In this paper, using these great advantages of GNR in the MIR wavelengths, a structure is proposed to be simultaneously used as a filter, demultiplexer, and power splitter, which could reduce the device area for some particular applications, where two of the devices may be needed simultaneously in associated working frequency. To this end, three different structures are presented based on GSNR resonators, numerically analyzed using the 3D-FDTD method. Depending on the number and location of output GNR waveguides (output port), these structures may have practical applications such as filter, wavelength demultiplexer, and power splitter, which all are fundamental elements in the nanophotonic and nanoplasmonic integrated circuits. In all of these structures, GNRs with a width of a few nanometers are used, which makes the structures presented simple in terms of design and fabrication, and most importantly, ultra-compact to be utilized in nanoscale devices at MIR functional region. The organization of the remains of the paper is as follows: in "The Structural Properties and Theoretical Analysis Method," the proposed structures and their various parameters are explained along with the theories and methods. "Results and Discussions" presents the results and discussions related to the proposed structures and finally, the conclusion is made in "Conclusion."

The Structural Properties and Theoretical Analysis Method
The schematic diagram of the proposed structures with the detailed geometrical parameters is shown in Fig. 1. All the three structures shown in this figure (a-c) use identical GNRs as incoming and outgoing waveguides, while a graphene-based square nanoring (GSNR) resonator completes the role of the filter (in Fig. 1a), demultiplexer (Fig. 1b), or the power splitter (Fig. 1c), depending on the location and the layout of the output GNR(s). The geometrical parameters include L 1 = 100 nm and W s = 20 nm (the side lengths and the width of the GSNR), L 2 = 100 nm and L 3 = 150 nm (the distances between the considered ports, that is where the detectors are placed, and the edge of the GNRs), W n = 20 nm (the width of the GNRs), and d = 5 nm (the coupling distances between the GNRs and the GSNR). The coupling distance, d, is chosen so that it guarantees an effective coupling of the surface plasmon polariton (SPP) waves between the components [41,42]. The width of the GNRs is fixed through all of the simulations to assure of the edge mode propagation only and the potential application of the structure in the ultra-compact integrated devices. The structures are shown in Fig. 1 are designed to be patterned on a sapphire film deposited on the silicon substrate. The Si/Al 2 O 3 substrate is chosen because it has strong absorption and low loss propagation in the MIR range. Therefore, it is the best candidate to support GNRs instead of SiO 2 . The refractive index and extinction coefficient k of Al 2 O 3 (Palik) are 1.59 and 0.0020 at 6 μm, respectively [43]. Also, in all the simulation runs, the thickness of the substrate is selected to be 50 nm for Al 2 O 3 and 500 nm for Si.
The design, dimensions, and materials of the proposed structures have been chosen to meet the fabrication considerations. The development of the proposed structure is technically feasible with well-developed patterning methods that have recently been used to produce GNRs with similar geometrical size features [44][45][46]. In a typical fabrication process, a high-quality large-area graphene layer can be grown on a copper substrate. After the etching process of Cu, it floats on the deionized water and, then, is pulled out on the substrate. Using electron-beam lithography and oxygen plasma, the GNRs can be patterned and the device is formed [47,48]. It is worth mentioning that for the proposed structures, the graphene monolayer is considered to be transferred on the Si/Al 2 O 3 substrate.
Graphene, as an ultra-thin film with a thickness of Δ ≈ 0.34 nm, guides the surface plasmonic waves in the MIR region. The surface conductivity of a graphene single layer is yielded from the well-known Kubo formula [34,49,50]: where ω, μ c , ℏ, T, k B , e, and τ are the angular frequency of the electromagnetic wave, chemical potential of the graphene, the reduced Plank's constant, working temperature, the Boltzmann's constant, the electron charge, and the carrier relaxation time of the material, respectively. Two terms are contained in the Kubo formula including the intra-band (the first term in Eq. 1) and the inter-band (the second term) transition contributions. In the MIR wavelengths, which is the objective working region in this article, the energy of photons, E ph , is low enough which let the μ c be always above half of it (i.e., ℏω < < 2μ c ), then the inter-band transition can be ignored here, and the σ g is reduced to the semi-classical Drude model [34,51]: determines the absolute value of the graphene chemical potential. The momentum relaxation time depends on the DC mobility (μ), the μ c , and the Fermi velocity (v F ≈ 10 6 m/s) of the material; Through the simulations, the single-layer graphene is modeled as an anisotropic dielectric, with the in-plane permittivity. If ε 0 is to be the vacuum permittivity, then the effective permittivity of the material is expressed as Eq. 4 [34,51].
The resonance wavelength of the GSNR, which works as the resonator through all the simulations, is [26,34]: Here, the integer m is the order of the resonance mode and L eff is the effective side length of the resonator, which is the mean of the internal and external sides of the GSNR resonator. In the analysis of the performance of the proposed structures, a key parameter is the n eff a demonstration of the effective refractive indices in the SPP modes which are supported by the GSNR. n eff is affected by several parameters within the models which are explained in the next section, and it can be obtained from the dispersion equation with the η 0 (≈377 Ω) as the impedance of the air [52]: Therefore, the propagation constant of the SPPs in the GSNR is β SPP = k 0 n eff , where k 0 is the free space wave vector.
The models are simulated using 3D-FDTD in Lumerical FDTD Software package with 24 identical perfectly matched layers as the absorbing boundary conditions all around the numerical simulation zone. The dipole point source, 2 nm over the location of the input monitor (P 1 ), was used to excite the incident SPP waves with transverse magnetic (TM) polarization and the power of P in , in the + x direction. The output monitors (P 2 -P 4 ) are set at the end of the output GNR waveguides to detect the transmitted powers P out . The normalized transmittance is then calculated as T = P out / P in . As long as edge modes are activated, the orientation of the dipole has an insignificant impact on the results. To obtain reliable, accurate, and at the same time, time-effective solutions, the numerical calculations used non-uniform mesh grid sizes of 0.1 nm for graphene, in the z-direction, while in the vicinity of the graphene components, the mesh size was 1 nm × 1 nm in the x-y plane.

Results and Discussions
In this section, the functions of the different structures presented in Fig. 1 are introduced and analyzed, such as the filter, demultiplexer, and power splitter. In all simulations, μ c for GNR waveguides is assumed to be constant and equal to 0.4 eV, unless otherwise is expressed. Since the output can be tuned using μ c adjustment, the 0.4 eV is chosen for this parameter as the mean value. Higher or lower chemical potentials are also applicable which may change the output spectrum/transmission. The μ c for the GSNR can be flexibly adjusted by controlling different gate voltages applied to graphene.

Filter
To analyze the performance of the plasmonic band-pass filter (BPF), the structure shown in Fig. 1a is proposed. This structure consists of a GSNR and two GNR waveguides as input/output waveguides. The dipole point source excites the SPP waves at the edge of the input GNR waveguide (edge mode), and they are then coupled to the GSNR resonator. These waves are backward and forward waves at the edge of the input GNR waveguide as well as the GSNR resonator, causing standing waves in the GSNR (here GSNR acts like the Fabry-Perot resonator (FPR)). As a result, the input wavelengths which satisfy the GSNR resonance conditions (Eq. 4) can be effectively coupled to the output GNR and the other wavelengths are stopped. In Fig. 1a, the two monitors are positioned P 1 and P 2 , respectively, to detect incident power and transmitted power, and the transmittance is defined as T = P 2 / P 1 . The transmission spectrum of the structure is shown in Fig. 2a. The geometrical parameters of the structure are equal to L 1 = L 2 = 100 nm, W s = W n = 20 nm, and d = 5 nm. The μ c for GNR waveguides and GSNR is considered constant at 0.4 eV. Three transmission peaks in the wavelength range of 6.5 to 11 μm pass through the structure, which corresponds to the resonant modes of the 4th order (λ 4 = 7.39 μm) and the 3rd order (λ 3 = 8.26 μm), and the 2nd order (λ 2 = 10.11 μm), respectively. These peaks have transmission amplitudes of 0.44, 0.54, and 0.58, respectively. To evaluate the performance of the BPF, we use the quality factor, which is defined as Q = λ 0 /Δλ, where λ 0 and Δλ are the central wavelength and full width at half maximum (FWHM), respectively. Table 1 shows the specifications of the different orders of resonance of the designed structure. As shown in Table 1, the quality factor decreases with an increasing central wavelength of the resonant mode. Therefore, there is a trade-off between these two parameters. Table 2 shows the comparison between the maximum transmission ratio in this work with other graphene-based structures, which compared to other structures, the proposed structure has a larger maximum transmission ratio. Figure 2 (b)-(d) show the spatial distribution profiles of H z from the cross-sectional view xy for the wavelengths corresponding to the fourth, third, and second-order resonant modes, respectively. The GSNR resonator also acts as an FPR cavity.
The effect of different parameters such as μ c , side lengths of GSNR (L 1 ), coupling distance between input/output GNR waveguides with GSNR (d), and width of GSNR (W s ) on the filter transmission spectrum is shown in Fig. 3.
The transmission spectra calculated for the different values of the μ c of the GSNR obtained from applying the external gate voltage are shown in Fig. 3a. The other parameters are set as shown in Fig. 2. As can be seen from Fig. 3a, the resonance wavelength of the GSNR decreases with increasing μ c , and the transmission spectrum experiences a blueshift. According to Eqs. (2, 5 and 6), this blueshift can be explained because increasing μ c leads to a reduction in n eff and thus a decrease in the resonance wavelength. Moreover, with increasing μ c , transmittance power increases. Figure 3b shows the effect of changing the L 1 on the transmission spectrum. It is clear that by increasing L 1 for a constant μ c (0.4 eV), the transmission peaks undergo redshift, which is resulted in Eq. 5. Besides, the coupling distance (d) between the input/output GNR waveguides and the GSNR affects the transmission peaks and the resonance wavelength, as shown in Fig. 3c. In this figure, by increasing the parameter d, besides the transmission spectrum experiencing a small blueshift resulting from the phase shift of the SPPs, the transmission peak and the stop-band level also decreases. This is due to the weakening of the coupling distance between the waveguide and the resonator. Furthermore, the transmission spectrum of the W s is investigated in Fig. 3d. As W s increases from 15 to 20 nm, the transmission spectrum shows a blueshift because the n eff of the GSNR decreases, and according to Eq. 4, the resonance wavelength decreases. It is clear from the simulation results that the structure presented in Fig. 1(a) can be used as a plasmonic BPF with dynamic tunability without re-fabrication of the structure in the MIR Fig. 2 (a) The transmission spectrum of the proposed BPF structure is shown in Fig. 1 spectral region. On the other hand, the proposed structure can be utilized in plasmonic devices due to its simplicity in terms of fabrication. Figure 1b shows the structure of a plasmonic demultiplexer consisting of a GSNR which is performed as a resonator and four GNR waveguides as input/output waveguides which are symmetrically placed around the GSNR. The geometrical parameters of the structure are similar to structure Fig. 1a, and since the proposed structure is symmetrical, only port 1 which is set as input, and ports 2, and 3 which are set as two outputs (output port 3 has the same position as port 4) are examined as follows, and the results are the same for the two output ports 3 and 4 due to symmetry. To reveal the incident power and transmitted power, three monitors are used in positions P 1 , P 2 , and P 3 , respectively, and the transmittance is defined as T = P 2 / P 1 and T = P 3 / P 1 , and the distance between these three monitors and the GSNR is equal to L 2 = 100 nm. The demultiplexer transmission spectra with the parameters μ c = 0.4 eV, L 1 = L 2 = 100 nm, W s = W n = 20 nm, and d = 5 nm are shown in Fig. 4a. In transmission spectra, graphene SPP waves can only be transmitted through the output port 2 during the 3rd resonance wave (λ = 8.29 μm) and have a high transmission amplitude of 0.56, while the transmission from through output port 3 (vertical port) is  . 3 The transmission spectra of the proposed BPF structure are shown in Fig. 1(a) for different a μc, b side lengths (L1) of GSNR, c coupling distance (d) between input/output GNR waveguides and GSNR, and d width (Ws) of GSNR negligible and equal to 0.01. This result completely agrees with the rule of the signal intersection at the resonance wavelength (port 1 ← → port 2 and port 3 ← → port 4). For the 2nd resonance wavelength (λ = 10.18 μm) and the 4th resonance wavelength (λ = 74.74 μm), the amplitudes of the transmission for the two output ports are approximately equal to 0.18 and 0.15, respectively. Hence, the structure is simultaneously served as two devices consisting of a demultiplexer (3rd resonance wavelength), and a tunable 1 × 3 coupler (2nd and 4th resonance wavelengths) with a constant splitting ratio (1:1:1) or 1 × 2 coupler when port 4 is closed. For better understanding, the transmission mechanism of the proposed structure, the spatial distribution profile of H z for the wavelengths corresponding to the 4th, 3rd, and 2nd order resonant modes are depicted in Fig. 4b-d, respectively. As shown in Fig. 4b, d, the edge modes are excited and propagated evenly along with the output of the GNR waveguides, and easily form a tunable 1 × 3 coupler. It is also clear from Fig. 4c that for output port 2, the resonant mode is located in the GSNR resonator, and the exciting edge mode and power pass-through output port 2, and almost no power can be coupled to output ports 3 and 4. Moreover, the demultiplexer performance of the structure can be evaluated by two important factors, namely, insertion loss (IL) and extinction ratio (ER). These two factors are defined as IL = 10 log(P w /P i ) and ER = 10 log(P w /P uw ), respectively [57], in which P i represents the power at the input waveguide, and P w and P uw are related to the desired and undesirable optical powers in the output waveguides. For the 3rd resonance wavelength (λ = 8.29 μm), the output port 2 and the output port 3 represent the desired (P w ) and undesirable (P uw ) output powers, respectively. Therefore, by applying the above equations, IL and ER for the 3rd resonance wavelength are equal to 2.5 dB and 16.55 dB, respectively.

Wavelength Demultiplexer
The performance of various parameters in this structure is shown in Fig. 5. The demultiplexer transmission spectrum for different μ c from 0.38 to 0.42 eV with 0.02 steps with geometrical parameters L 1 = L 2 = 100 nm, W s = W n = 20 nm, and d = 5 nm is illustrated in Fig. 5a. It is obvious in Fig. 5a that with an increase in μ c , the transmission spectrum Fig. 4 (a) The transmission spectra of ports 2 and 3 of the proposed wavelength demultiplexer structure are shown in Fig. 1(b). (b)-(d) The spatial distributions of Hz at the peaks of the transmission for resonance modes experiences a blueshift, as explained in the previous section. The effect of GSNR dimensions on the transmission spectrum and its results are shown in Fig. 5b. Other parameters are constant and equal to μ c = 0.4 eV, W s = W n = 20 nm, and d = 5 nm. As it is known, by increasing the dimensions of the GSNR, the transmission spectrum undergoes a redshift, which can be explained according to the GSNR conditions (Eq. (4)); since the GSNR dimension increases, the resonance wavelength increases. Another important parameter effect on the transmission spectrum of the two output ports 2 and 3 is the coupling distance between the input/ output GNR waveguides and the GSNR, which is shown in Fig. 5c. The μ c for GSNR and GNR waveguides is the same and is equal to 0.4 eV. As can be seen from this figure, the coupling length between the GNR waveguides and the GSNR has a low effect on the resonance wavelength; because the resonant mode changes only with the geometrical parameters and μ c of the GSNR (transmission spectrum has a negligible blueshift due to SPP phase shift). It is also clear that as the coupling length increases, the coupling loss increases as well, thus the transmission amplitude drops. To have ultra-compact devices with a high ER, the width of the GNR input/output waveguides and the GSNR must be chosen so that the dimensions of the device are compressed and reduced. Therefore, in all structures, the width of the GNR and GSNR is less than the critical values to be able to excite only the edge mode. However, the effect of different W s on the transmission spectrum of the two output ports 2 and 3 is shown in Fig. 5d. By increasing the W s from 15 to 25 nm, the transmission spectrum of the two output ports experiences a redshift, in addition to increasing the amplitude. This behavior can be understated by using Eq. (4) and is the same as the previous section.

Power Splitter
To employ the effective coupling between the output GNR waveguides and the GSNR, two GNR waveguides are added  Fig. 1(b) for different a μc, b side lengths (L1) of GSNR, c coupling distance (d) between input/output GNR waveguides and GSNR, and d width (Ws) of GSNR in parallel at the top and bottom of the GSNR, as shown in Fig. 1c. By adding these two GNR waveguides to the output, the input energy can have two practical waveguides. This structure can be used as a power splitter at some resonance wavelengths. The power splitter is one of the basic functional elements in nanoplasmonic integrated devices and circuits. In Fig. 1c, one monitor is placed in position P 1 to detect incident power and two other monitors are located on positions P 2 and P 3 to detect transmitted power, and the transmittance is defined as T = P 2 / P 1 and T = P 3 / P 1 . Here, only the effects of two parameters, namely, μ c and coupling distance (d) between two output GNR waveguides with GSNR are investigated, and the results are given as follows. The other parameters are fixed, similar to the previous two structures.
By applying the gate voltage and its effect on the μ c at the two output ports, the power splitting ratio can be altered, and thus the transmission spectrum changes at the two output ports, as is shown in Fig. 6a. Here, μ c of the output port 2 is fixed to μ c = 0.40 eV, and the μ c of the output port 3 is set once to 0.25 eV and, then, to μ c = 0.10 eV. When the μ c of output ports 2 and 3 is , the transmission amplitude for output port 2 is 0.42, and for the output port 3 is 0.12. The spatial distribution profile of H z for this case is shown in Fig. 6d. In the 3rd resonance wavelength (λ = 8.29 μm), the transmission for the output port 3 is 0.34 and for the output port 2 is 0.25, and the spatial distribution profile of H z is given in Fig. 6c. In these two resonance wavelengths, the structure acts as a power splitter that has different splitting ratios. These results can be justified using the graphene conductivity equation Eq. 2 because increasing the μ c leads to increasing the σ g . Therefore, the transmission loss along with the GNR decreases and the transmission amplitude at the output ports increases with increasing μ c . In contrast to the above discussion, when the wavelength is 8.29 μm, the transmission at output port 3 is slightly more than the transmission at output port 2. Since the coupling between output waveguide 2 and GSNR is slightly weaker than output port 3. Therefore, the input signal beam is coupled directly to the output port 3. For the 4th resonance wavelength (λ = 7.67 μm), the value of the transmission amplitude for the output port 2 is 0.31, while the value for the output port 3 sharply drops to near zero (0.02). Therefore, the proposed structure acts as a demultiplexer. The spatial distribution profile of H z for this case is shown in Fig. 6b. When the μ c of the output port 3 is adjusted to 0.10 eV, the SPP transmission loss in the output port 3 is greatly increased, the SPP waves are reflected, and it cannot support any graphene surface plasmon in the wavelength range of 6 to 12 μm.
The spatial profile of H z in Fig. 6e confirms this issue. The spatial profile of H z for μ c2 = μ c3 = 0.40 eV, is shown in Fig. 6f, and the transmission amplitude for both ports are the same and equal to 0.42.
The transmission spectrum of the power splitter structure in Fig. 1c strongly depends on the coupling distance (d) between the two output GNR waveguides and the GSNR. The strong optical coupling between the output GNR waveguides, and the GSNR requires a small mismatch in the propagation constant of the SPPs along with the GNR, which satisfies the phase-matching conditions. The coupling ability between the output GNR waveguides and the GSNR can be managed by changing the d, as shown in Fig. 7a. In this figure, the transmission spectrum is shown for a situation where the coupling distance between the GSNR and output port 2 and port 3 is 5 nm and 10 nm, respectively. For the 3rd wavelength (λ = 8.31 μm), the transmission for both output ports is 0.21, and the proposed structure operates similar to a 1 × 2 power splitter. The spatial distribution profile of H z for this case is shown in Fig. 6c. For the 2nd resonance wavelengths (λ = 10.11 μm) and the 4th resonance wavelengths (λ = 7.44 μm), the structure acts as a demultiplexer. At these resonance wavelengths, the transmission amplitude for output port 2 is 0.44 and 0.33, respectively, while the transmission amplitude for output port 3 is strongly suppressed and is nearly zero. The spatial distribution profile of H z for these two resonance wavelengths is given in Fig. 6b, d. This figure shows that as the length of the coupling increases, the intensity of the coupling to the output port decreases dramatically.

Conclusion
In summary, three fundamental devices including a filter, a demultiplexer, and a power splitter were designed for applications in nanophotonic integrated circuits. The abovementioned arrangements were theoretically and numerically analyzed using the 3D-FDTD method. These structures consisted of the main block, i.e., a GSNR, as a resonator, and a set of some extra GNR waveguides as input/outputs. The structures were very convenient and simple in terms of implementation and practical applications. The coupling intensity of GSPPs between GSNR and GNR output waveguides is affected by various parameters such as chemical potential, coupling length, and dimensions of GSNR, and so does the transmission spectrum at the output of waveguides, as a consequent. Changing these parameters, the coupling intensities could be controlled and tuned in the MIR spectral range.
Author Contribution Design, methodology, and numerical simulations of the paper: Morteza Janfaza. Writing of the manuscript: Morteza Janfaza. Writing-comments and suggestions: Mohammad Ali Mansouri-Birjandi and Alireza Tavousi.

Availability of Data and Material
All datasets generated and/or analyzed during this study are available from the corresponding author on reasonable request.

Code Availability
The code used during the current study is available from the corresponding author on reasonable request.

Declarations
Ethics Approval The authors have followed the ethical principles and accurate references to scientific sources in the original paper.