The ground-state properties of silicene are calculated using standard density functional theory (DFT) computations implemented in the Abinit software [23], specifically within the local density approximation (LDA) [24]. The Kohn-Sham (KS) electron wave functions are expanded in the plane-wave ( P W ) basis [25], represented as: P W k + G r = Ω 0 − 1 / 2 e i k + G ∙ r (1) where k is a wave vector in the first Brillouin zone along the Γ Κ Μ Γ path, G is a reciprocal lattice vector, and Ω 0 represents the unit-cell volume associated with the real-space lattice. The number of P W s is limited by the energy cut-off k + G / 2 ≤ 25 Hartree ( ∼ 680 eV). To remove the effect of core electrons, norm-conserving pseudopotentials of the Troullier-Martins type are utilized [26]. The 3D periodic boundary conditions required for plane-wave DFT calculations are achieved by replicating the silicene sheets along the z direction with a vacuum distance of 20 Å . The scope of this study is limited to ideal freestanding silicene, and geometric optimization was unnecessary as its impact on the estimation of charge carrier velocity was found to be negligible beyond the third significant digit (Supplementary Figure S1). The lattice constant and buckling are fixed at 3.82 Å and 0.45 Å , respectively.

To calculate the electron band structure of silicene, we prepared two datasets of parameters. The first set involves a 90 × 90 × 1 Monkhorst-Pack grid [27] that is Γ-centered and unshifted, which can yield energy levels up to 10 eV above the Fermi level using 50 bands. The second set uses a high-resolution grid of 540 × 540 × 1 , which enables the calculation of the linear band structure of silicene near the K point and includes up to 6 bands due to computational constraints. The KS electronic structure, which includes the ground-state energy and electronic density, is then utilized as the starting point in the GW calculation.

To obtain an accurate calculation of the band dispersion of freestanding silicene, it is essential to incorporate many-body effects, which can be achieved by using the many-body GW self-energy method. As well-known, the GW method is a widely used approach for improving the accuracy of DFT calculations. The self-energy in the GW approximation is given by the expression [28]: ∑ r , r ′ , ω = i 2 π ∫ − ∞ ∞ d ω ′   G   r , r ′ , ω − ω ′ W r , r ′ , ω ′ (2) where G is the one-particle Green’s function and W is the screened Coulomb interaction. Eq. 2 represents the product of Green’s function G and the dynamically screened interaction W ω . The screened Coulomb interaction is given by: W ω = ϵ − 1 ω ∙ υ (3) where ϵ − 1 ω is the inverse of the dynamical dielectric function, and υ is the bare Coulomb interaction. To make our analysis simpler, we choose to ignore vertex corrections in both the self-energy and the polarizability. Instead, we construct the G and W components in the GW method using the DFT-KS electronic structure. To integrate over frequency, we use the contour deformation (CD) method [29], which is a reliable and precise method of frequency integration. This involves deforming the real axis contour to calculate the self-energy as an integral along the imaginary axis while accounting for a contribution from the residual of the contour-included poles of G . This method is essential for our analysis as it ensures the accuracy of our calculations.

To explore the plasmon characteristics of SiNSs that have a width of 100 nm or greater, we adopt the approach introduced by [21]. The investigation involves determining the plasmon frequency ( ω p ), which is accomplished through the utilization of the following expression: ω p = R e 2   π   e 2   N 2 D ϵ   m ∗ q cos 2 ⁡ θ − v 2 4 − i v 2 (4) the parameters of Eq. 4 are detailed as follows.

• e is the conventional electron charge

The charge density ( N 2 D ) can also be expressed as a function of the one-dimensional (1D) charge density ( N 1 D ) of an isolated nanostrip: N 2 D = N 1 D d (5) where d is the vacuum distance between contiguous strips.

From the experimental results presented by [30], it has been demonstrated that in graphene nanostrips with widths ranging from 155 to 480 nm, there exist two distinct resonance modes: the surface plasmon and the edge plasmon. It is also expected that these modes are present in SiNSs. In particular, the edge plasmon mode can be selectively tuned by altering the ribbon width [31], while the surface plasmon mode is found to be more sensitive to various external factors, including doping levels, temperature, electron mobility, and the angle of plasmon excitation [32]. These findings have significant implications for the design of nanoscale electronic and photonic devices. Indeed, by controlling the dimensions of the nanostrips and selectively tuning the plasmon modes, it may be possible to engineer novel devices with tailored optical and electronic properties.

Eq. 4 is a mathematical expression that describes only the dispersion of the surface plasmon mode. While this equation does not account for the nature of the edge plasmon mode, it is still an efficient means of calculating the frequency dispersion of surface plasmons. Furthermore, this expression is consistent with experimental observations, which have shown that the plasmon wavelength follows the sample length, with the sample length being much larger than both the vacuum distance and the ribbon width [33]. For simplicity, we use the term “plasmon” exclusively to refer to the surface plasmon mode in this work.

To further customize the analysis of plasmonic properties in a specific context, Eq. 4 can be modified as needed, for instance, using.

(i) The Fermi level ( E F ) shift [34] as:

In Eq. 7, it is assumed that E F is in the middle of the bandgap, where E c , E v , N c and N v are the conduction band edge, the valence band edge, the effective DOS in the conduction band, and the effective DOS in the valence bands, respectively. The other well-known parameters k B and T represents the Boltzmann constant and the absolute temperature. Thus, Eq. 4 can be expressed in terms of the needed Fermi level to inject or eject electrons as well as the desired temperature.

Due to Eq. 4 being a straightforward analytical expression, when the plasmon damping ( γ ) is very high, the radicand becomes negative. The plasmon damping is identified [36] as: γ = ν 2 (8)

From the physical point of view, this effect can be caused by various mechanisms, including scattering, absorption, and radiation, which lead to a loss of plasmon energy. As a result, the plasmon response can move towards larger momenta, which corresponds to a higher frequency or shorter wavelength.

Moreover, the notion of a complex dielectric function ( ϵ M ) [37] can be used to better comprehend this phenomenon, which is represented by: ϵ M q , ω = ϵ R ω + i ϵ I ω (9) thus, ϵ R ω corresponds to the real part of the in-phase response, while i ϵ I ω corresponds to the imaginary part of the out-of-phase response. Eq. 9 indicates that the plasmons can be observed by detecting zeros in the real part within a frequency range where the imaginary part is minimum or cero. Indeed, forbidden plasmon regions appear if this condition is not met.

With this in mind, it is now feasible to compute the effective electron mass in Eq. 4 by incorporating the charge-carrier velocity ( v F ) of silicene, as outlined in [21]: m ∗ = ∆ 2   v F 2 (10) and the bandgap ( ∆ ) of a single nanostrip can be evaluated as [21]: ∆ = 2   π   v F   ℏ w (11) here, ℏ is the Planck constant.

As noted, the charge-carrier velocity of silicene is the critical parameter in Eq. 4 and is the foundation of the semi-analytical model. Then, to estimate the charge-carrier velocity accurately, we employed either GW or LDA-DFT calculations to perform a linear fit of the π (highest valence) and π * (lowest conduction) bands. This linear fit is equivalent to the Dirac cone approximation [38], denoted as: E = ± v F p (12) here, ± sign represents the conduction and valence bands and p is the momentum represented as p = ℏ k − Γ Κ .

On the other hand, to estimate the band structure of SiNSs, it is crucial to take into account the quasi-one-dimensional confinement of charge carriers. This confinement leads to the formation of multiple energy sub-bands ( E n ), expressed as [21, 39]: E n = ± ∆ 2 n 2 + 2 p ∥ 2 m ∗ ∆ (13) here the integer number n represents the sub-band index ( n = 1 , 2 , 3 , … ) and p ‖ is the parallel wavevector to the nanostrip direction. To further clarify, in quasi-1D systems such as SiNSs, the confinement of charge carriers leads to quantization of the energy levels, resulting in the formation of multiple sub-bands. These sub-bands are distinguished by their energy levels and wavefunctions and can be indexed by n . More importantly, Eq. 13 displays a parabolic band dispersion at the Γ point, which in fact is predominant in narrow-wide SiNSs in the eV scale, and starts to be negligible by increasing the ribbon width ( w → ∞ ), i.e., the unique electronic properties of silicene are recovered gradually.

Now, to show the plasmon spectrum (i.e., the maximum of the plasmon peak) for selected values of q , the standard approach of spectral line profile can be used: L = 1 1 + 4 ω − ω o 2 W 2 (14) here, L is the Lorentzian function, ω o is the transition energy/frequency of the maximum (taken from data using Eq. 4), ω is the energy/frequency sampling of interest, and capital W is the full width at half maximum (FWHM). FWHM is set to 0.25 for all spectra.

The Lorentzian line function is a widely used model for describing spectroscopic features in physical systems like ions, atoms, and molecules, as well as in SiNSs. This study employs frequency units for ω , ω o and W , following the customary practice for collective excitations spectra [40]. It is worth noting that in samples supported on metallic substrates, the maximum of the plasmon peak can shift due to core-electron excitations [41], but this is not the case for freestanding SiNSs. Furthermore, the plasmon structure in SiNSs is not expected to be a simple Lorentzian peak; however, this approach effectively demonstrates the control of plasmons in SiNSs.

To begin, we analyze the band structure of freestanding silicene using a 90 × 90 × 1 grid. Figure 1 shows a comparison between the band structure obtained from LDA-DFT (dashed black line) and quasiparticle GW (solid red line) calculations within an energy range of ± 4 eV. Although the impact of exchange and correlation effects is evident, we observed only minor changes in the band shapes and band dispersions. Interestingly, regardless of the theoretical level used, we found that silicene displays a conical band dispersion at the Κ point, which is commonly referred to as a Dirac cone.

It is worth noting that the LDA results show a lowering of the bands, while the gaps at the Γ (3.3 → 4.1 eV) and M (1.7 → 2.21 eV) points are increased by up to 19% and 23%, respectively. Such behavior is a typical effect of GW approximation, as observed in previous studies. Moreover, there is an increase in the linear energy dispersion around the K point (GW results), which confirms a rise in the charge-carrier velocity of freestanding silicene. It is interesting to note that this phenomenon has also been observed in freestanding graphene when comparing the two theoretical levels [29].

We now focus on the region around the Dirac cone (K point) by examining energy-momentum data obtained from GW calculations of freestanding silicene. We have considered the first conduction ( π * ) band and last valence ( π ) band with the second (high resolution) grid of 540 × 540 × 1 , as shown in Figure 2A. By applying a linear fit (cyan line) to this data, we can estimate the charge-carrier velocity. Interestingly, we observe a perfect linear dispersion of electron energy within a certain range of energy ( ± 0.16 eV) and wave vector ( ± 0.04 Å − 1 ), which closely matches the Dirac cone approximation (as described in Eq. 12).

Figure 2B displays the linear behavior of charge-carrier velocity, although only within a narrower energy range ( ± 0.13 eV) since the intense (white) lines start to disappear. These findings are significant as they provide limits for the semi-analytical model used (∼31 THz), which relies on a linear dependence. Similar results are obtained through LDA-DFT calculations (results not shown here), suggesting a minor limit of approximately ± 0.11 eV (∼26 THz). This information is crucial for understanding the plasmonic properties of SiNSs at the THz frequencies.

Specifically, Table 1 compares the estimated values by LDA-DFT ( v F = 0.540 × 10 6 m s⁻¹) and GW ( v F = 0.742 × 10 6 m s⁻¹) calculations. The LDA-DFT calculation underestimated the GW value by 27.22%. Currently, there is no consensus on the experimental charge-carrier velocity of silicene. To gain a deeper understanding of our findings, we compared the results obtained for silicene with those of graphene at the same theoretical levels (LDA-DFT-graphene (Supplementary Figure S2) and GW-graphene (Supplementary Figure S3). Interestingly, we found that the charge-carrier velocity of silicene is 33.63% lower than that of graphene. It is worth noting that the GW calculation for graphene showed only a slight overestimation of 1.64%, which closely matched the experimental result for graphene (Exp-graphene) [42]. Similarly, as in silicene, LDA-DFT calculations for graphene yield a charge-carrier velocity that is 24.64% lower than the experimentally measured value. This discrepancy highlights the limitations of the LDA-DFT method and the need to employ more advanced calculations. Additionally, these results confirm that the GW approximation is a reliable method for predicting the charge-carrier velocity of silicene, and future experiments are likely to produce similar outcomes.

Equation 11 demonstrates that there is an inverse relationship between the bandgap and ribbon width. Specifically, as the width of the ribbon increases, the bandgap value decreases exponentially. Therefore, if the ribbon width were to approach infinity ( w → ∞ ), the bandgap value would approach zero ( Δ → 0 ), similar to silicene. This implies that the structural and electronic properties of silicene would be gradually restored. Figure 3 illustrates this relationship between ribbon width and bandgap, demonstrating the exponential decrease in bandgap value, which remains consistent regardless of whether the charge-carrier velocity was determined through GW (red curve) or LDA-DFT (blue curve).

The impact of ribbon width on the bandgap is most pronounced apparently in narrow SiNSs with widths ranging from 10 to 50 nm, as shown in Figure 3A, where the bandgap decreases from 0.3 eV to 50 meV. Conversely, wider SiNSs with widths ranging from 100 to 500 nm, shown in Figure 2B, exhibit a decrease in bandgap from 30 meV to 5 meV. In both cases, the bandgap decreases by a factor of six. Moving forward, we will focus on the electronic properties of SiNSs using the charge-carrier velocity of silicene obtained through GW calculations (i.e., v F = 0.702 × 10 6 m s^-1).

Table 2 presents information on the bandgaps and effective electron masses for various nanostrips with different widths ranging from 100 to 500 nm. The data reveals a noteworthy trend: the wider the nanostrip, the smaller the bandgap. This trend is evident in the 100 nm wide nanostrip, which has a bandgap of approximately 30 meV, compared to the 500 nm wide nanostrip, which has a bandgap of approximately 6 meV. Furthermore, when comparing the same nanostrips, the effective electron masses show the same trend. The effective electron mass values obtained for the different nanostrips are consistent with both experimental findings and predictions for graphene nanostrips [43], [45], in terms of orders of magnitude. These results suggest that the properties of the SiNSs are similar to those of graphene nanostrips.

In Figure 4, we now present the band structure (Figures 4A, C) and density of states (DOS) (Figures 4B, D) of the 100 nm and 500 nm wide nanostrips, in a ± 70 meV energy range. This analysis helps in understanding the electronic properties of these nanostrips and their potential applications in various fields. As an example, SiNSs with a bandgap of 1.1 eV, similar to silicon, could be useful for the semiconductor industry. On the other hand, nanostrips with a bandgap in the range of a few meV, particularly SiNSs with widths equal to or greater than 100 nm wide (see Table 2), are highly promising for nanophotonics applications.

Upon comparison, it becomes apparent that there is a quadratic band dispersion of the conduction and valence bands near the Γ point (zero) for both 100- and 500-nm wide systems. It should be noted that we have maintained the same effective electron mass for all bands, with only the sub-band index, n , being altered in Eq. 10. Though it is possible to accurately calculate the effective electron mass of each band, including both conduction and valence bands, using ab initio methods such as DFT, the main issue is that the systems analyzed in this study are too broad for current computers to handle. Then, the semi-analytical model proves advantageous in this regard and produces reasonable outcomes. In light of this, we can observe that in the same energy range, the 500 nm wide nanostrip (Figure 4C) displays a considerably greater number of bands in comparison to the 100 nm wide nanostrip (Figure 4D). It is also worth noting that with an increase in the width of nanostrips, the energy gap between the bands decreases, and they tend to converge, whenever possible, to the same energy state, similar to what we observe in silicene with the corresponding π * and π bands (Figure 1, Figure 2A).

Figures 4B, D show the DOS of the 100 nm and 500 nm wide nanostrips, respectively. These figures display several peaks corresponding to the bands in the band structure plot. Interestingly, smoothing the DOS histogram (red line) can lead to a narrower bandgap of the system. For example, the bandgap of the 100 nm wide nanostrip is roughly 23 meV, while that of the 500 nm wide nanostrip is around 4 meV. This results in a 23% and 33% reduction in the bandgap, respectively (Table 2).

In this section, we investigate the plasmonic properties of SiNSs, focusing on the effect of charge carrier density ( N 2 D ). Charge density in a nanostrip refers to the amount of electric charge per unit area. There are several methods to estimate charge density experimentally or theoretically, including analyzing the diffraction pattern of X-rays or neutrons scattered by a crystal, using electrostatic force microscopy, performing Hall effect measurements, and measuring the capacitance of a semiconductor device. DFT calculations can also be used to determine charge density.

Indeed, previous studies have reported that the charge carrier density in isolated graphene nanostrips is approximately N 1 D ≈ 10 15 cm⁻¹ [39],[46]. Similar results are expected for SiNSs. However, this value can vary depending on factors such as doping, ribbon size, or geometric factors. Then, Eq. 5 can be used as a starting point to adjust the charge density to the desired value (see Table 3). Interestingly, as the distance between adjacent nanostrips increases, the charge density decreases. For instance, when the vacuum distance between nanostrips increases from 10 nm to 100 nm, the charge density reduces by one order of magnitude. However, it is worth noting that, in this case, the Fermi level ( E F ) only decreases by approximately threefold, going from 0.38 eV (for d = 10 nm) to 0.12 eV (for d = 100 nm).

Figure 5A shows the dispersion of plasmon frequency-momentum for a 100 nm wide nanostrip as a function of the reciprocal wave vector ( q ). As expected for 2D materials such as graphene and silicene, the plasmon dispersion in SiNSs follows a q -like behavior, which can be explained by the nature of Eq. 4. Interestingly, increasing the distance ( d ) between adjacent nanostrips from 10 nm to 100 nm leads to a decrease in the plasmon frequency and plasmon dispersion. This phenomenon occurs because the plasmon frequency is proportional to the square root of the charge density ( ω p ∼ N 2 D ).

Additionally, Table 3 indicates that a separation distance of 10 nm results in a high charge density ( N 2 D = 1 × 10 13 cm⁻²) and a shift in the Fermi level of 0.38 eV, which may lead to excessive doping for a 100 nm wide nanostrip with a small bandgap of a few meV (∼30 meV). Therefore, we focus on a charge density of N 2 D = 1 × 10 12 cm⁻², which is a commonly used value for modeling graphene or similar materials. This choice results in a Fermi level shift of only 0.12 eV and a separation distance of 100 nm is adequate. Furthermore, this vacuum distance ensures that the system can be treated as a quasi-isolated SiNS.

Considering this, Figure 5B shows the dispersion of the maximum of the plasmon peak for a 100 nm wide nanostrip with a charge density of 1 × 10 12 cm⁻² at selected momentum ( q ) values within a frequency range of 0.1–3 THz, where in turn, the most relevant silicene-based plasmonic applications, could occur. On the other hand, Figure 5C presents a complete analysis of the effect of the charge density from 1 × 10 12 to 1 × 10 13 cm⁻². The analysis reveals two significant observations: i) forbidden regions for plasmons at values of q near zero (indicated by the purple region), and ii) the spectral weight enhances with increasing charge density and momentum. At high charge densities ( ∼ 10 13 ), a plasmon response can be detected up to approximately 35 THz (red region).

The effective electron mass ( m * ) is a key parameter being studied. It characterizes the behavior of electrons in a solid and is dependent on the bandgap and charge-carrier velocity. Experimental methods such as transport measurements, the Hall effect, and angle-resolved photoemission spectroscopy (ARPES) can be used to determine the effective mass. Alternatively, quantum mechanical models such as DFT, tight-binding models, and effective mass approximation (EMA) can be used for its calculation.

Eq. 10 contains this parameter, which is only affected by changes in the bandgap, such as those resulting from modifications to the ribbon width (Table 2) since the charge-carrier velocity ( v F = 0.702 × 10 6 m s⁻¹) remains constant in this work. With this in mind, Figure 6A shows the plasmon dispersion of SiNSs with widths ranging from 100 to 500 nm as a function of the momentum ( q ). Increasing the width of the nanostrip causes the effective electron mass to decrease, increasing the plasmon frequency. This correlation is straightforward since the plasmon frequency is inversely proportional to the square root of the effective electron mass, i.e., ω p ∼ 1 / m * . This trend is consistent with the values in Table 2, where the effective electron mass decreases as the nanostrip width increases.

Figure 6B provides compelling evidence of the controllability and tunability of plasmon response, showing that increasing the width of the nanostrip or decreasing the effective electron mass results in an increase in plasmon frequency. This figure focuses on a frequency range of 0.1–3 THz, which is highly relevant for plasmonic applications in silicene, as mentioned earlier. Therefore, nanostrip systems ranging from 100 to 500 nm hold promise as potential candidates for such applications. Furthermore, Figure 6C demonstrates that there are no forbidden regions for plasmons and confirms that a decrease in the effective electron mass leads to an increase in plasmon frequency by approximately 25 THz (red region).

A final remark, Fei et al. [30] have successfully prepared graphene nanostrips with widths similar to those examined in this study (i.e., w = 155, 270, 380, and 480 nm).

The plasmon excitation angle ( θ ) is determined by the interaction of electromagnetic waves with a material and depends on factors such as energy, momentum, material properties, and system geometry. ARPES can experimentally measure the plasmon excitation angle by directly measuring the momentum and energy of electrons in the material. Hence, Figure 7A explores the impact of the angle of plasmon excitation on a 100 nm wide nanostrip, which can be altered in the experimental setup. It is worth noting that as the angle of plasmon excitation increases, the frequency of the plasmon decreases, indicating that at an angle of 90°, the plasmonic response to the perpendicular direction becomes non-existent. It is evident that for angles less than 90°, the plasmonic response can be customized for specific purposes within the frequency range of 0.1–3 THz (Figure 7B) by adjusting this parameter. In addition, increasing the angle shifts the plasmon response towards lower frequency values.

Figure 7C exhibits the existence of forbidden regions where plasmons are prohibited, and as the angle increases. Interestingly, these prohibited regions can also be observed by increasing the momentum value (as evidenced by the extended purple region). Moreover, it is important to highlight that the highest spectral intensity is achieved with small angles ( ω p ∼ cos 2 θ ), producing a plasmonic response of nearly 10 THz (red region).

Electron relaxation rate ( ν ) is the rate at which electrons in a material lose their energy. Time-resolved pump-probe spectroscopy and transient absorption spectroscopy are experimental techniques that can measure the relaxation rate of electrons in a material. DFT calculations or many-body perturbation theory can also be used to estimate this parameter.

Then, we focus on the electron relaxation rate in nanostrip of 100 and 500 nm wide. Additionally, this parameter is crucial in providing the values for electronic mobility and plasmon relaxation rate (Eq. 8) as shown in Table 4. Particularly, the electron mobility ( υ ) can be obtained using the relationship υ = e   τ / m * , where τ represents the electron relaxation time, which includes the Fermi level and other parameters mentioned in Eqs 6, 7.

Figure 8A depicts the plasmon frequency-momentum dispersion of a 100 nm wide nanostrip as a function of the increasing electron relaxation rate. Notably, the curves appear to overlap, prompting us to investigate whether this effect is attributable to the width of the strip. To that end, we analyzed a 500 nm wide nanostrip (as shown in Figure 8B) but found that the plasmon frequency-momentum dispersion followed the same pattern.

Given that this effect is not due to strip width, we continue the analysis to a lower momentum regime for a 100 nm wide nanostrip, as depicted in Figure 9. We observe in Figure 9A that the plasmon frequency continues to be proportional to the square root of q . However, when the relaxation rate of electrons increases, we observe two additional phenomena.

(i) The plasmonic response moves towards higher momentum values.

Apart from the significant influence of the electron relaxation rate, Figure 9B shows that even a 100 nm wide SiNS subjected to the electron relaxation rate effect displays plasmon responses in the frequency range of 0.1–3 THz, suggesting even greater adjustable plasmonic properties. Figure 9C demonstrates the existence of forbidden regions for plasmons, which become more prominent as the electron relaxation rate and momentum increase (purple region). As previously mentioned, increasing the electron relaxation rate results in a significant reduction in the plasmon frequency. For instance, in a 100 nm wide nanostrip, the plasmon frequency reaches a maximum of around 2 THz (red region).

Hence, Eq. 4 provides a comprehensive understanding of the parameters that influence the plasmon response, including its frequency and dispersion. By thoroughly examining these parameters, one can manipulate them individually or in combination to customize the plasmon response for a specific application.