Quenching of low-energy optical absorption in bilayer C₃N polytypes

Many recent experimental and theoretical works have focused on a detailed understanding of the electronic and optical properties of two-dimensional (2D) materials, due to their potential use in the design of innovative optoelectronic devices which could combine atomic-size dimensions with improved performance Mak et al. (2010); Zhang (2015); Chernikov et al. (2014); Cudazzo et al. (2016); Lin et al. (2016); Rousseau et al. (2021); Zhang et al. (2017); Galvani et al. (2016); Tran et al. (2014); Thygesen (2017); Palummo et al. (2015). Great attention has been reserved to the study of Van der Waals structures, where atomically-thin 2D materials are vertically stacked and held together by weak and long range dispersive interactions, strongly affecting both the optical Andersen et al. (2015); Paleari et al. (2018) and the electronic properties of the isolated single layers He et al. (2014). A promising material for such novel applications is graphene-like 2D polyaniline (also known as monolayer C₃N), which has been recently synthesized through different approaches Mahmood et al. (2016); Yang et al. (2017): its electronic and optical properties have been extensively studied from a theoretical point of view Wu et al. (2018); Bonacci et al. (2022); Tang et al. (2022); Zanfrognini et al. (2023), revealing a quasi-particle band structure with indirect band-gap and intense optical absorption for photons in a narrow spectral range around 2 eV.

Vertical stacking of two layers of C₃N is a possible approach to tune its electronic and optical properties. Few theoretical works have analysed the stability of bilayers C₃N (BL-C₃N) as a function of the possible stacking patterns Wu et al. (2018); Wenya et al. (2021). Following the notation of Ref. [Wenya et al., 2021], all these calculations have obtained negative formation energies for AB and AB^′ (displaced-like), as well as AA^′ (sandwich-like stacking) arrangements. Furthermore, BL-C₃N with AA^′ and AB^′ stackings have also been experimentally synthesized, as described in Ref. [Wenya et al., 2021], where a detailed investigation of the electronic properties using scanning tunnelling spectroscopy (STS) has been provided. The results indicate a strong change of the electronic-transport band gap passing from monolayer to bilayer, together with relevant modifications of the electronic properties as a function of the stacking sequence.

Motivated by these experimental advances, in this work we discuss the optical properties of BL-C₃N via first principles methods, properly including excitonic effects which are known to play a fundamental role in the description of optical absorption in 2D-materials Latini et al. (2015); Andersen et al. (2015); Hüser et al. (2013). As a first result, we find for all systems a strong optical absorption in an energy region around 1.7 eV which is sightly lower than the absorption peak of the isolated monolayer Bonacci et al. (2022), despite the consistent reduction of the electronic gap. This feature, combined with the energy band structure associated to these absorption peaks, is promising for photovoltaics. Moreover, our results indicate that, for all the considered stackings, the optical response does not present other bright excitonic states at lower energies. Such behaviour is not observed in other common semiconducting bilayer homo-structures, e.g. BL-hBNPaleari et al. (2018), BL-MoS₂Mak et al. (2010); Palummo et al. (2015) or BL-phosphoreneTran et al. (2014). Finally, the origin of this low energy absorption quenching is explained and rationalized, focusing on the properties of the single particle states involved in the formation of the lowest-energy excitons.

The article is organized as follows. In Section II, we summarize the computational methods used within this work, while in Section III we present quasi-particle band structures for the AB and AB^′ displaced-like stacking patterns. In Section IV we discuss the optical properties of these bilayers and in Section V we provide a rationale for the negligible oscillator strength observed for the low-lying excitons. Finally, in Section VI, we discuss the optical properties of bilayer C₃N with AA^′ sandwich-like stacking, predicting a strong optical absorption quenching starting from the symmetry properties of this polytype and confirming such interpretation with an approximate solution of the Bethe Salpeter equation.

Ground state structural and electronic properties have been investigated using density functional theory (DFT), as implemented in the plane-wave-basis-set package Quantum ESPRESSO Giannozzi et al. (2009, 2017). In these calculations, we have used norm-conserving ONCV pseudopotentials Hamann (2013), within the GGA-PBE approximation Perdew et al. (1996) for the exchange-correlation potential. Van der Waals interactions between layers have been taken into account by adding the dispersion correction proposed by Grimme Grimme (2006) to the exchange-correlation energy (PBE-D2). Equilibrium structural properties have been obtained by relaxing both the in-plane unit cell and the atomic positions up to when the components of the forces acting on each atom were smaller than 5$\cdot$10^-4 Ry/Bohr. In all ground state calculations we used a 12$\times$12$\times$1 Monkhorst-Pack Monkhorst and Pack (1976) k-grid to sample the Brillouin Zone (BZ) and a kinetic energy cutoff of 90 Ry for the plane wave basis set used to represent single particle wavefunctions.

In the case of AB and AB^′ stacking motifs, Kohn-Sham wavefunctions and eigenvalues, computed from the equilibrium ground state charge density, have been used to evaluate quasi-particle (QP) corrections to DFT energies within the GW approximation Hedin (1965); Onida et al. (2002) for the electron self energy. QP corrections $\epsilon_{n\mathbf{k}}^{\text{KS}}$ have been computed using the single shot G₀W₀ approach, evaluating the expectation value of the operator ${\Sigma}-{v}_{\text{KS}}$ on the Kohn-Sham states $|\psi_{n\mathbf{k}}\rangle$, being ${\Sigma}$ the electron self-energy operator and ${v}_{KS}$ the exchange-correlation potential. When solding the Dyson equation we linearized the frequency dependence of the self-energy, ${\Sigma}(E^{\text{QP}}_{n\mathbf{k}})\approx{\Sigma}(\epsilon^{\text{KS}}_{n\mathbf{k}})+\frac{\partial{\Sigma}}{\partial\epsilon}\big{|}_{\epsilon^{\text{KS}}_{n\mathbf{k}}}(E^{\text{QP}}_{n\mathbf{k}}-\epsilon^{\text{KS}}_{n\mathbf{k}})$, as proposed in Ref. [Hybertsen and Louie, 1986].

The electron-electron screened interaction $W$ has been computed using the Random Phase Approximation (RPA), as implemented in the Yambo code Marini et al. (2009); Sangalli et al. (2019). Converged QP gaps within a 10 meV threshold required the inclusion of 700 bands and a $\mathbf{G}$-vector cutoff of 16 Ry in the construction of the screening matrix at the RPA level. The frequency dependence of $W$ has been described using the Godby-Needs plasmon-pole model Godby and Needs (1989), and 1000 bands have been included in the sum-over-states appearing in the correlation part of the self energy. To reduce spurious interactions among different cells in the stacking direction, we have used a supercell length along $z$ of 23.5 Å, together with a 2D cutoff Ismail-Beigi (2006); Rozzi et al. (2006) for the Coulomb potential. Finally, to speed-up the convergence of QP gaps w.r.t. the k-point mesh, we have adopted the approach recently proposed by Guandalini et al. Guandalini et al. (2022). In this work we have verified that, with this method, a 18$\times$18$\times$1 Monkhorst-Pack k-grid already provides converged gaps within the chosen threshold of 10 meV.

Starting from QP corrected electronic energies, we have obtained excitonic properties by solving the Bethe-Salpeter equation (BSE) Strinati (1988); Rohlfing and Louie (2000) in the resonant (Tamm-Dancoff) approximation, i.e. via diagonalization of the excitonic Hamiltonian

where $K^{d}$ ($K^{x}$) is the direct (exchange) part of the BSE kernel, $(v,v^{\prime})$ and $(c,c^{\prime})$ are the valence and conduction bands included in the BSE and $E_{(v,c)\mathbf{k}}^{QP}$ are the QP corrected electron energies. Converged exciton energies have been obtained including the two highest-occupied valence and the four lowest unoccupied conduction bands in the construction of $H_{exc}\big{(}vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k^{\prime}}\big{)}$, and using a 48$\times$48$\times$1 Monkhorst-Pack k-grid to sample the BZ. We point out that, in the solution of the BSE, QP corrections have been approximated via a scissor-stretching operator (see Supplemental Material sup for a detailed description about the fitting procedure), while the electron-electron screened interaction has been computed at the RPA level, in the static approximation, using the same converged parameters adopted for the calculation of QP corrections.

Finally, by diagonalization of $H_{exc}\big{(}vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k^{\prime}}\big{)}$

we obtained the energies $E_{\lambda}$ and the envelope functions $A_{\lambda}\big{(}vc\mathbf{k}\big{)}$ for each exciton $\lambda$. Optical absorption is finally computed as the imaginary part of the macroscopic dielectric function,

In the above expression, $V$ is the unit cell volume and $D_{\lambda}$ the oscillator strength of exciton $\lambda$ defined as

where $\hat{\epsilon}$ is the in-plane polarization direction of the incoming light, and $d_{vc\mathbf{k}}=\langle\varphi_{v\mathbf{k}}|\mathbf{r}|\varphi_{c\mathbf{k}}\rangle$ the single-particle interband dipole matrix element.

In Figure 1 we present the crystal structures of BL-C₃N with AB (upper panel) and AB^′ (lower panel) stackings: yellow (light blue) spheres denote Carbon (Nitrogen) atoms and small (large) atoms are located on the upper (lower) layer, denoted from now on as $\mathrm{L}_{1}$ ($\mathrm{L}_{2}$). For both stacking motifs, we have obtained an in-plane lattice parameter of 4.849 Å, slightly smaller than that of the isolated monolayer (4.857 Å); the interlayer distance (evaluated as the separation along $z$ between Carbon atoms with the same in-plane coordinates) has been found equal to 3.22 Å for the AB and 3.21 Å  for the AB^′ stacking. These values are in agreement with those obtained with PBE-D2 calculations in Ref. [Wenya et al., 2021], while slightly smaller than the interlayer distances computed with VdW-functionals in Ref. [Wu et al., 2018].

We now briefly discuss the crystal symmetries of the two stackings. The point group of AB-C₃N is D_3d, and also includes non-symmorphic operations. This stacking possesses a spatial inversion center (red dot in the upper panel of Fig. 1) together with a three-fold C₃ rotation axis along $z$-direction, whose in-plane position is denoted by a green dot in Fig. 1. Furthermore, this stacking motif is invariant under mirror reflections $\sigma$ w.r.t. planes parallel to the stacking direction and represented by dashed red lines in Fig. 1: These planes are respectively denoted as $\hat{\sigma}_{\Gamma\mathrm{M}}$, $\hat{\sigma}_{\Gamma\mathrm{M^{\prime}}}$, and $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$ as they are aligned to these high symmetry directions in the BZ. AB^′-C₃N has lower symmetry: Its point group is C_2h, that includes the inversion symmetry, a two-fold in-plane rotation axis lying between the two C₃N planes (represented by the green dashed line in the lower panel of Fig. 1 and denoted as $\hat{\mathrm{C}}_{2}^{\Gamma\mathrm{K}}$), and a mirror symmetry plane $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$ parallel to the $\Gamma\mathrm{M^{\prime\prime}}$ direction in the BZ.

In Figure 2 we present the electronic band structure of AB-C₃N (panel a) and AB^′-C₃N (panel b), computed first within DFT-PBE (dashed green lines) and then including QP corrections at the G₀W₀ level (solid black lines). We note that the two stackings are characterized by an indirect band gap, both at the DFT and GW level. In the case of AB stacking, the electronic band dispersion along $\Gamma$M and $\Gamma\mathrm{M^{\prime\prime}}$ coincide, as a consequence of the three-fold rotational symmetry. Therefore, the highest-energy valence band has six equivalent maxima, while the bottom of the conduction band is located at the $\Gamma$ point. We also note that the presence of doubly degenerate bands at $\Gamma$ is consistent with irreducible representations of dimension 2 in the $\mathrm{D}_{3\mathrm{d}}$ point group. The calculated indirect gap at PBE level is 0.108 eV, which is increased to 0.72 eV once QP corrections are taken into account. Finally, the direct band gap is found along the $\Gamma\mathrm{M}$ direction, at k located approximately at half the $\Gamma\rightarrow\mathrm{M}$ path: the G₀W₀ gap is 1.85 eV, 0.73 eV larger than the value obtained at the PBE level (1.12 eV).

In the AB^′ stacking, the maximum of the valence band is found at the $\mathrm{M}^{\prime\prime}$ point, while the lowest unoccupied conduction state is at $\Gamma$, as in the case of the AB stacking motif. Our calculations give an indirect gap of 0.136 eV at the DFT-PBE level, which is increased to 0.73 eV with the inclusion of QP corrections. The obtained indirect gap for AB^′-C₃N is slightly smaller than the one measured experimentally with STS on SiO₂-Si substrates in Ref. [Wenya et al., 2021] (0.85$\pm$0.03 eV), but in better agreement than other GW calculations Wu et al. (2018) where the Hybertsen-Louie plasmon-pole model was used Hybertsen and Louie (1986). Because of the lack of three-fold rotational symmetry around the $z$ axis, the directions $\Gamma\mathrm{M}$ and $\Gamma\mathrm{M}^{\prime\prime}$ are no longer equivalent. As a result, the minimum direct gap is found approximately at half the $\Gamma\rightarrow\mathrm{M}^{\prime\prime}$ path, with a value around 1.79 eV (1.09 eV within DFT) while the indirect gap between the conduction at $\Gamma$ and the top-valence at $\mathrm{M}$ is slightly larger than the one between $\Gamma$ and $\mathrm{M}^{\prime\prime}$. We obtain a $\Gamma_{c}-\mathrm{M}_{v}$ gap of 0.87 eV (0.249 eV) at the G₀W₀ (DFT-PBE) level. Electronic gaps computed within the G₀W₀ approximation for both stacking motifs at high symmetry points are summarized in Table 1.

Comparing the band dispersions for AB and AB^′ stackings, we notice that, for $\mathbf{k}$ along the $\Gamma\mathrm{M}$ direction, the lowest pair of conduction bands are almost degenerate (splitting of about 1 meV) in the AB^′ stacking, while well separated in the AB case (splitting larger than 0.2 eV). In the Supplemental Material sup , we provide a qualitative explanation of this peculiar lack of splitting among the two lowest conduction bands in AB^′-C₃N, analysing the quasi-symmetries of the sublattice where conduction states are localized once $\mathbf{k}$ is taken along $\Gamma\mathrm{M}$ direction. A similar behaviour is also observed for the two topmost valence states, where the splitting is, however, not negligible also in the AB^′ stacking (splitting of about 50 meV). Overall, the obtained gaps for the AB and AB^′ stacking are similar (in agreement with the hybrid-DFT results of Ref. [Wenya et al., 2021]), while the electronic dispersions differ because of the different symmetry properties of the two stackings. We point out that there are clear differences with the monolayer case, as we observe a strong reduction in both the indirect and direct electronic gaps (see Table. 1 for a comparison).

We complete our analysis by recalling that both stacking configurations are invariant under mirror reflection $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$, therefore electronic states for k along this direction can be classified in terms of their parity w.r.t. such symmetry. We have numerically found that, in both stackings, the two highest occupied valence bands are even w.r.t. $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$, (here denoted as $v_{i}^{+}$ in Fig. 2, with $i=1,2$). On the other hand, if k has modulus in the range $[\frac{|\Gamma\mathrm{M}|}{3},|\Gamma\mathrm{M}|]$ the two lowest unoccupied conduction bands are $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$-odd ($c_{i}^{-}$). Instead, for k close to $\Gamma$, the second and the third conduction bands are even w.r.t $\hat{\sigma}_{\Gamma\mathrm{M^{\prime\prime}}}$ and exhibit a stronger dispersion with k, if compared with odd conduction states. We point out that a similar analysis can be carried out also for electronic states along this direction in monolayer C₃N. In that case, the highest valence (lowest conduction) is even (odd) w.r.t. the mirror symmetry along $\Gamma\mathrm{M}^{\prime\prime}$ direction, while the second conduction is even. In the case of AB stacking, the same parity analysis can be presented for the bands along the direction $\Gamma\mathrm{M}$, as the crystal is invariant w.r.t. $\hat{\sigma}_{\Gamma\mathrm{M}}$ mirror symmetry. Such symmetry classification will be exploited to understand bilayer optical properties in the following.

We now turn to the discussion of the optical properties of bilayer C₃N. In Figure 3 we show the absorption spectra computed for AB (panel a) and AB^′ bilayer C₃N (panel b), at the independent particle level (dotted black lines), and with the inclusion of electron-hole interaction by solving the BSE as detailed in Section II. Green dashed (solid red) lines have been obtained assuming light polarized along the ${\Gamma\mathrm{K}}$ (${\Gamma\mathrm{M}^{{}^{\prime\prime}}}$) direction. For clarity, polarization versors $\hat{\epsilon}$ are shown together with the crystal structures in the insets.

The AB spectrum is dominated by an intense peak (denoted as C) at energy $E\approx 1.70$ eV, whose spectral position and intensity are not dependent on the polarization direction. Such peak is due to a set of almost degenerate eigenstates of H_exc, characterized by relevant contributions from single particle transitions between the valence band $v^{+}_{1}$ and the conduction state $c^{-}_{1}$ along $\Gamma\mathrm{M}$ equivalent directions. Among these, transitions with the highest weights are denoted by arrows on the band structure shown in Fig. 3(c). We point out that transitions between $v^{+}_{2}$ and $c^{-}_{2}$ for $\mathbf{k}$ along the same direction also contribute to this absorption peak, though with a smaller weight than $v^{+}_{1}\rightarrow c^{-}_{1}$, and are omitted for clarity in Fig. 3(c). The situation is slightly different in the case of AB^′ C₃N. Indeed, also for this stacking motif the absorption spectrum is dominated by a single intense peak, but its position in energy and its strength depend on the in-plane light polarization direction. In Fig. 3(b), we have labelled as C₁ the main peak at 1.71 eV obtained for light polarization along ${\Gamma\mathrm{K}}$ and as C₂ the absorption maximum at 1.73 eV found for ${\Gamma\mathrm{M}^{{}^{\prime\prime}}}$-polarized light. The C₁ peak is mainly due to $v^{+}_{1}\rightarrow c^{-}_{1}$ for $\mathbf{k}$ along $\Gamma\mathrm{M}^{{}^{\prime\prime}}$ direction, as shown schematically by green arrows in Fig. 3(d), again with a smaller contribution coming from $v^{+}_{2}$ to $c^{-}_{2}$ transitions for the same $\mathbf{k}$ points (not shown in Fig. 3d). On the other hand, the C₂ absorption peak comes from transitions between the two highest occupied valence states and the two lowest (quasi-degenerate) conduction bands along ${\Gamma\mathrm{M}}$ and ${\Gamma\mathrm{M}}^{\prime}$ directions (see red arrows in Fig. 3d).

First, we note that the energy window found for strong optical absorption in both stackings is still in a good range for solar energy conversion, not much lower than that found for monolayer C₃N of 1.82 eV in Ref. [Bonacci et al., 2022]. Furthermore, the associated transitions are still in region favourable to electron-hole splitting, as in the monolayer. This can be particularly interesting, since one can expect the bilayer to have a more stable structure, compared to a monolayer, when deposited on a suitable electrode.

Second, and central to this work, the striking feature of bilayer C₃N optical spectra is the apparent absence of intense absorption peaks due to strongly bound excitonic states formed by single particle transitions close to the electronic direct gap. As shown in the insets of Fig. 3(a,b), within the energy range between 1.25 and 1.5 eV, two absorption structures (labelled in both cases as P₁ and P₂) are present, but they exhibit optical strengths which are almost two orders of magnitude smaller than the most intense peaks.

In AB-C₃N, both P₁ and P₂ peaks are due to a pair of degenerate excitons, respectively at energies $E_{\mathrm{P}_{1}}$ = 1.35 eV and $E_{\mathrm{P}_{2}}$ = 1.47 eV, with oscillator strengths not dependent on the polarization direction. For completeness, we point out that diagonalization of the excitonic Hamiltonian also provides other excitonic resonances (the lowest with energy of 1.34 eV) which are characterized by null oscillator strength within numerical accuracy. Such excitons, dark by strict symmetry reasons, will not be considered further in the following. P₁ and P₂ excitons are almost totally composed by electron-hole transitions between the highest occupied valence $v^{+}_{2}$ and the lowest unoccupied conduction $c^{-}_{1}$ states, with wave-vectors $\mathbf{k}$ along $\Gamma\mathrm{M}$ and equivalent directions in the BZ. This is better clarified in Fig. 4(a,c), where we show plots of the quantity

for the P₁ and P₂ excitons, respectively. In Eq. (5), the index $v$ ($c$) is fixed to the last valence (first conduction) band and the sum over $\lambda$ is performed over the pair of degenerate states responsible for the P₁ and P₂ peaks. We notice that single particle transitions forming the excitons P₁ are mainly localized in the central region of $\Gamma\mathrm{M}$ and equivalent directions, with $A^{vc}(\mathbf{k})$ having non-negligible values for $|\mathbf{k}|$ mainly in the interval $[\frac{1}{3},\frac{2}{3}]|\Gamma\mathrm{M}|$. On the other hand, excitons P₂ (Fig. 4c) are still localized along $\Gamma\mathrm{M}$ directions, but the corresponding function $A^{vc}(\mathbf{k})$ has intense contributions from transitions slightly closer to the $\mathrm{M}$ point and exhibits a node for $\mathbf{k}$ points along this direction.

In the case of the AB^′ stacking, P₁ and P₂ peaks are related each to a single-nondegenerate exciton at energies $E_{\mathrm{P}_{1}}$ = 1.30 eV and $E_{\mathrm{P}_{2}}$ = 1.42 eV, respectively. Differently from the AB case, the oscillator strength of these excitations is polarization dependent as shown in the inset of Fig. 3(b), where we notice that both resonances are dark for light polarization along the ${\Gamma\mathrm{M}}^{\prime\prime}$ direction while they exhibit a small but non-zero optical activity for incoming light polarized along ${\Gamma\mathrm{K}}$. As in the AB stacking, such excitons are mainly formed by transitions between the highest occupied valence and the lowest conduction band. In Fig. 4(b,d) we show the functions $A^{vc}(\mathbf{k})$ computed for exciton P₁ and P₂, respectively. We see that, in both cases, the transitions involved in these excitations are strongly localized along the single $\Gamma\mathrm{M}^{\prime\prime}$ direction of the BZ (see e.g. Fig. 2), with the P₁ exciton having the main contributions coming from the middle of $\Gamma\mathrm{M}^{\prime\prime}$, where the minimum direct electronic band gap is located, and the P₂ resonance characterized by a node along this BZ-line.

The observed polarization dependence in the optical absorption of AB^′-C₃N can be rationalized via symmetry arguments, similarly to the analysis proposed in Ref. [Paleari et al., 2018] for bilayer hBN. The point group of AB^′-C₃N is C_2h, so the in-plane exciton dipole operator projected along ${\Gamma\mathrm{M}}^{\prime\prime}$ (${\mathbf{D}}_{\Gamma\mathrm{M}^{\prime\prime}}$) transforms as the B_u irreducible representation of C_2h, while ${\mathbf{D}}_{\Gamma\mathrm{K}}$ as A_u. In fact, ${\mathbf{D}}_{\Gamma\mathrm{M}^{\prime\prime}}$ is invariant under the mirror symmetry $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$ (as oriented along $\Gamma\mathrm{M}^{\prime\prime}$), but it changes sign under inversion and $\hat{C}^{\Gamma\mathrm{K}}_{2}$ rotation (see Fig. 1). Differently, ${\mathbf{D}}_{\Gamma\mathrm{K}}$ is even w.r.t. $\hat{C}^{\Gamma\mathrm{K}}_{2}$ and odd under $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$ and inversion, therefore behaving as the $A_{u}$(C_2h) irreducible representation. Therefore, the eigenstates $|\lambda\rangle$ of $H_{exc}$ such that $\langle 0|\mathbf{D}|\lambda\rangle\neq 0$ ($|0\rangle$ being the excitonic vacuum transforming as the fully symmetric representation A_g) transform as A_u or B_u if ${\mathbf{D}}$ is projected along ${\Gamma\mathrm{K}}$ or ${\Gamma\mathrm{M}}^{\prime\prime}$, respectively. Considering for example A_u states, these will have null optical activity by symmetry, once incoming light is polarized along ${\Gamma\mathrm{M}}^{\prime\prime}$. Therefore, we can explain the presence of absorption peaks in AB^′-C₃N which turn on and off according to the chosen polarization direction. Furthermore, following the presented analysis, we can assign P₁ and P₂ excitons to the A_u representation, exactly as the eigenstates of H_exc responsible for the C₁ peak, while C₂ is expected to be due to excitations transforming as B_u.

The situation is different for AB stacked C₃N. For simplicity, we discuss the brightness of excitonic eigenstates using the subgroup C_3v of D_3d, formed by the mirror planes $\sigma$ along $\Gamma\mathrm{M}$ directions together with the two 3-fold rotations depicted in the upper panel of Fig. 1. The in-plane exciton dipole operator $\mathbf{D}$ transforms as the two-dimensional irreducible representation E of C_3v. As a consequence, all bright excitons behave as E of C_3v, so that they are expected to be double-degenerate and characterized by an isotropic oscillator strength, in agreement with the numerical results obtained from ab initio calculations.

We now turn our attention to the quasi-dark nature of the low-lying bound excitons P₁ and P₂ in both the considered stackings. We point out that such small optical activity cannot be related to an interlayer nature of these excitons, i.e. it is not due to a negligible spatial overlap between electron and hole wavefunctions. This point is better clarified by looking, for example, at the real space wavefunctions of exciton P₁, as shown in Fig. 5, with AB (AB^′) results reported in the left (right) panel. For each stacking, the upper (lower) wavefunction has been obtained assuming the hole (represented by the red dot) fixed on layer L₁ (L₂) and located in the plane close to a nitrogen atom.

The wavefunctions clearly indicate the intralayer nature of such excitons: In fact, looking at the excitonic wavefunction isosurfaces, we notice that the electron has a high probability to be found on the same layer on which the hole is localized.

Therefore, the negligible dipole strength of these low-lying excitons is a consequence of the small interband dipole associated to the electron-hole single particle transitions involved in these excitations. This can be understood from the independent particle (IP) absorption spectrum shown in Fig. 3(a,b) as dotted lines, where, in both cases, the optical signal is negligible for energies close to the direct QP gap, with the IP absorption onset located at higher photon energies. In the following section we discuss in greater detail the single particle states involved in these low-lying excitons and we propose a possible rationale for the observed IP absorption quenching.

In this Section we develop a model for the electronic bands in proximity of the direct gap to rationalize the small oscillator strength of the single particle transitions close to the direct electronic gap. As the following analysis is valid for both stackings, here we focus on the case of AB^′-C₃N (results for AB motif are presented in the Supplemental Material sup ). We will restrict our analysis to $\mathbf{k}$ along $\Gamma\mathrm{M}^{\prime\prime}$, in the region where valence to conduction transitions giving the highest contribution to the low-energy quasi-dark excitons are located.

As already discussed in the literature Wenya et al. (2021); Wu et al. (2018), the lowest lying conduction bands and the highest valence states have a $\pi$ character. therefore, we can analyse them using a tight-binding (TB) Hamiltonian, obtained considering one 2p_z orbital for each atom. In the following we will denote as $\mathbf{\tau}_{\alpha}$ the position of both the atom $\alpha$ and the 2p_z orbital localized on it. In practice, we construct a TB Hamiltonian as

where $\mathbf{R}$ is a lattice vector and $t(\alpha 0,\beta\mathbf{R})=\langle\alpha 0|{H}|\beta\mathbf{R}\rangle$ are the hopping matrix elements between two 2p_z orbitals, localized at sites $\mathbf{\tau}_{\alpha}$ and $\mathbf{\tau}_{\beta}+\mathbf{R}$, respectively. These matrix elements have been computed fully ab initio by Wannierization Souza et al. (2001); Marzari et al. (2012) of DFT bands using the Wannier90 Mostofi et al. (2008, 2014); Pizzi et al. (2020) code. Details about the procedure are provided in the Supplemental Material sup .

As there are 16 atoms in the unit cell, at a general $\mathbf{k}$ the TB Hamiltonian can be written as a 16$\times$16 Hermitian matrix, in a block-like form as

where $\mathbf{H}^{\mathrm{L}_{1}}$, $\mathbf{H}^{\mathrm{L}_{2}}$, and $\mathbf{H}^{\mathrm{IL}}$ are 8$\times$8 matrices corresponding to the different layers and their coupling. Indeed, to obtain this separation, we have grouped together the orbitals localized on L₁ and on L₂, associating to each orbital $\alpha_{1}$ localized at $\mathbf{\tau}_{\alpha_{1}}$ on layer L₁ an orbital $\alpha_{2}$ at $\mathbf{\tau}_{\alpha_{2}}$ on layer L₂, such that $\mathbf{\tau}_{\alpha_{2}}=\hat{\mathrm{I}}\mathbf{\tau}_{\alpha_{2}}$, $\hat{\mathrm{I}}$ being the inversion symmetry operator. Then $\mathbf{H}^{\mathrm{L}_{1}}$ ($\mathbf{H}^{\mathrm{L}_{2}}$) is the block of $\mathbf{H}^{\mathrm{2L}}$ which contains the intralayer hopping within layer L₁ (L₂), while $\mathbf{H}^{\mathrm{IL}}$ depends on the hopping integrals between orbitals on different layers. We can now write $\mathbf{H}^{\mathrm{2L}}(\mathbf{k})=\mathbf{h}_{\mathrm{IN}}(\mathbf{k})+\mathbf{h}_{\mathrm{IL}}(\mathbf{k})$, where $\mathbf{h}_{\mathrm{IN}}(\mathbf{k})$ is block-diagonal while $\mathbf{h}_{\mathrm{IL}}(\mathbf{k})$ is purely off-diagonal. We note that the subscripts ”IN” and ”IL” stand for intralayer and interlayer, respectively. Using spatial ($\hat{\mathrm{I}}$) and time ($\hat{\mathrm{T}}$) inversion symmetries, in Appendix A we show that, for each $\mathbf{k}$, $\mathbf{h}_{\mathrm{IN}}(\mathbf{k})$ has a spectrum of eigenvalues $\varepsilon_{n\mathbf{k}}^{0}$ which are degenerate in pairs. Furthermore, we also show that the eigenspace associated to each eigenvalue $\varepsilon_{n\mathbf{k}}^{0}$ is spanned by two Bloch states $|\phi^{\mathrm{L}_{1}}_{n\mathbf{k}}\rangle$ and $|\phi^{\mathrm{L}_{2}}_{n\mathbf{k}}\rangle$ such that $|\phi^{\mathrm{L}_{2}}_{n\mathbf{k}}\rangle=-\hat{\mathrm{I}}\cdot\hat{\mathrm{T}}|\phi^{\mathrm{L}_{1}}_{n\mathbf{k}}\rangle$ and with $|\phi^{\mathrm{L}_{i}}_{n\mathbf{k}}\rangle$ localized on layer $\mathrm{L}_{i}$.

Before proceeding, we clarify the physical meaning of the splitting of $\mathbf{H}^{\mathrm{2L}}$ into $\mathbf{h}_{\mathrm{IN}}$ and $\mathbf{h}_{\mathrm{IL}}$. In particular, $\mathbf{h}_{\mathrm{IN}}$ can be thought as an intralayer Hamiltonian, describing the two layers as not interacting with each other. Note however that the presence of the other layer is implicitly considered as the intralayer matrix elements in $\mathbf{H}^{\mathrm{L}_{1}}$ ($\mathbf{H}^{\mathrm{L}_{2}}$) are affected by the presence of opposite layer L₂ (L₁). On the other hand, $\mathbf{h}_{\mathrm{IL}}$ describes the interlayer interaction and acts as a perturbation of $\mathbf{h}_{\mathrm{IN}}$ since interlayer hopping integrals are typically smaller than the intralayer ones. With this interpretation, the states $|\phi^{\mathrm{L}_{1}}_{n\mathbf{k}}\rangle$ and $|\phi^{\mathrm{L}_{2}}_{n\mathbf{k}}\rangle$ can be thought as Bloch states with the same energy $\varepsilon_{n\mathbf{k}}^{0}$ localized on one of the two monolayers, if the coupling $\mathbf{h}_{IL}$ is set to zero. We now define $\varepsilon_{v\mathbf{k}}^{0}$ ($\varepsilon_{c\mathbf{k}}^{0}$) as the energy of the highest occupied valence (lowest unoccupied conduction) band on these two non-interacting layers. The effect of the interlayer coupling will be to mix these layer-localized wavefunctions, in order to give the electronic states of the bilayer.

While the discussion presented so far is general and valid for each $\mathbf{k}$-point of the BZ, we now specialize to $\mathbf{k}$-vectors along the $\Gamma\mathrm{M}^{\prime\prime}$ direction. For these wavevectors, $\mathbf{H}^{\mathrm{2L}}$ commutes with $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$ and the same is valid for $\mathbf{h}_{\mathrm{IN}}$ and $\mathbf{h}_{\mathrm{IL}}$ separately. Therefore,

with $i=1,2$ (we have also numerically verified these relations by computing the eigenstates of $\mathbf{h}_{\mathrm{IN}}$). We now include the effect of $\mathbf{h}_{\mathrm{IL}}$ using first order degenerate perturbation theory, separately diagonalizing the matrix representation of ${h}_{\mathrm{IL}}$ on the two subspaces $\{|\phi^{\mathrm{L}_{1}}_{v\mathbf{k}}\rangle,|\phi^{\mathrm{L}_{2}}_{v\mathbf{k}}\rangle\}$ and $\{|\phi^{\mathrm{L}_{1}}_{c\mathbf{k}}\rangle,|\phi^{\mathrm{L}_{2}}_{c\mathbf{k}}\rangle\}$. With this procedure, we obtain the two highest-energy ($\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$-even) valence bands $\{|\varphi_{v_{1}\mathbf{k}}\rangle,|\varphi_{v_{2}\mathbf{k}}\rangle\}$ and the two lowest-energy ($\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$-odd) conduction states $\{|\varphi_{c_{1}\mathbf{k}}\rangle,|\varphi_{c_{2}\mathbf{k}}\rangle\}$ in the bilayer. Such states have been labelled as ($v^{+}_{1}$,$v^{+}_{2}$) and ($c^{-}_{1}$,$c^{-}_{2}$), respectively, in Fig. 2(b). They can be compactly written as

In Eqs. (9)-(10), $n=v,c$, $j=1,2$, $s_{j}=-1$ ($+1$) for $j=1$ ($j=2$), $\Delta_{n\mathbf{k}}=\langle\phi^{\mathrm{L}_{1}}_{n\mathbf{k}}|{h}_{\mathrm{IL}}|\phi^{\mathrm{L}_{2}}_{n\mathbf{k}}\rangle$, and

where $\gamma_{n\mathbf{k}}=\mathrm{Arg}[\Delta_{n\mathbf{k}}]$ guarantees that the relative phase between the projections $c_{\alpha}(n_{j}\mathbf{k})=\langle\alpha\mathbf{k}|\varphi_{n_{j}\mathbf{k}}\rangle$ of a Bloch state $|\varphi_{n_{j}\mathbf{k}}\rangle$ on 2p_z orbitals localized on different layers is gauge invariant, i.e. it does not change under the transformation $c_{\alpha_{1}}(n\mathbf{k})\rightarrow e^{i\eta}c_{\alpha_{1}}(n\mathbf{k})$, $\eta$ being an arbitrary phase. See App. A for further details.

We point out that the zero-order expression given by Eq. (9) is a good approximation for the two lowest $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$-odd conduction bands and the two highest-energy $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$-even valence bands, for the considered $\mathbf{k}$ vectors along $\Gamma\mathrm{M}^{\prime\prime}$ direction. In principle one should also consider other terms in the expression of the eigenstatesSakurai and Napolitano (2017), coming from higher orders of the perturbative series, describing the coupling between $|\phi^{\mathrm{L}_{i}}_{n\mathbf{k}}\rangle$ and the eigenstates of $\mathbf{h}_{\mathrm{IN}}$ with different eigenvalues. For $\mathbf{k}$ points around the middle of the $\Gamma\mathrm{M}^{\prime\prime}$ direction, these terms can be neglected in first approximation, as the other eigenstates of $\mathbf{h}_{\mathrm{IN}}$ with the same $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime\prime}}$-parity of $|\phi^{\mathrm{L}_{i}}_{n\mathbf{k}}\rangle$ have energies far from $\varepsilon_{n\mathbf{k}}^{0}$ (w.r.t. to the interlayer coupling strength), so that the hybridization is negligible. Our numerical results (not shown) indicate that these neglected terms become more relevant for $\mathbf{k}$ along the same direction, but closer to the $\Gamma$ point. Therefore, looking at Eq. (9), we can understand that the lowest odd-conduction band and the highest even-valence can be seen, respectively, as antibonding and bonding combinations of the conduction and the valence states localized on the two monolayers. We recall that the states defined in Eq. (11) are still eigenstates of the intralayer Hamiltonian $\hat{h}_{\mathrm{IN}}$, with $|\tilde{\phi}^{\mathrm{L}_{2}}_{n\mathbf{k}}\rangle=-\hat{\mathrm{I}}\cdot\hat{\mathrm{T}}|\tilde{\phi}^{\mathrm{L}_{1}}_{n\mathbf{k}}\rangle$.

Starting from Eq. (9), we can evaluate the interband matrix element between the last occupied valence $|\varphi_{v_{2}\mathbf{k}}\rangle$ and the lowest unoccupied conduction $|\varphi_{c_{1}\mathbf{k}}\rangle$ as

where $\hat{\epsilon}$ is the light polarization. To make the treatment simpler, in the following we neglect the last two terms in Eq. (12), as they involve states localized on different layers, so that their value is generally small as a result of the reduced overlap among the wavefunctions. In this way, we obtain that the interband dipole is the difference between intralayer-interband dipoles:

The approach adopted to compute these matrix elements starting from the proposed tight binding model is discussed in the Supplemental Material sup .

If $\hat{\epsilon}$ is chosen along the $\Gamma\mathrm{M}^{\prime\prime}$ direction, we immediately find that $d_{\mathrm{L}_{1}}^{\epsilon}(\mathbf{k})=d_{\mathrm{L}_{2}}^{\epsilon}(\mathbf{k})=0$, because of the parity of layer-resolved states, given by Eq. (8). However, if $\hat{\epsilon}$ is taken along the direction $\Gamma\mathrm{K}$, orthogonal to $\Gamma\mathrm{M}^{\prime\prime}$, we cannot explain the quenching by straightforward symmetry arguments. To clarify this point, in Fig. 6(a) we show the modulus (upper panel) of the intralayer-interband dipoles $d_{\mathrm{L}_{1}}^{\Gamma\mathrm{K}}(\mathbf{k})$ and $d_{\mathrm{L}_{2}}^{\Gamma\mathrm{K}}(\mathbf{k})$, together with their relative phase (lower panel) $\Delta\varphi_{d}(\mathbf{k})=\mathrm{Arg}[d_{\mathrm{L}_{1}}^{\Gamma\mathrm{K}}(\mathbf{k})]-\mathrm{Arg}[d_{\mathrm{L}_{2}}^{\Gamma\mathrm{K}}(\mathbf{k})]$. Our results indicate that the obtained intralayer dipole matrix elements are equal in modulus and they exhibit a relative phase close to zero in the range of $\mathbf{k}$ points here considered, i.e. $|\mathbf{k}|$ in $[\frac{|\Gamma\mathrm{M}^{\prime\prime}|}{3},|\Gamma\mathrm{M}^{\prime\prime}|]$. As the total interband dipole is the difference among intralayer contributions – see Eq. (12) –, it will be almost zero, as a consequence of the destructive interference of the two layer-resolved components. In other words, the transition probability from $|\varphi_{v_{2}\mathbf{k}}\rangle$ to $|\varphi_{c_{1}\mathbf{k}}\rangle$ due to $\Gamma\mathrm{K}$-polarized light can be interpreted as the quantum superposition of the interband scattering processes occurring on the two layers separately, whose probability amplitudes are out-of-phase, giving an overall negligible interband oscillator strength.

In Fig. 6(b), the continuous red line indicates the interband dipole between $|\varphi_{v_{2}\mathbf{k}}\rangle$ and $|\varphi_{c_{1}\mathbf{k}}\rangle$, computed using the perturbative solution of the TB model (see Supplemental Material sup for details), while the red dots are the same quantities obtained ab initio using Yambo, to check the validity of our approximate treatment. We notice that, as this cancellation is not symmetry-constrained, the interband dipole is small, but not exactly zero. Such cancellation is exact at $\mathrm{M}^{\prime\prime}$ point, because of symmetry reasons. In fact, as $\mathrm{M}^{\prime\prime}$ is invariant under spatial inversion $\hat{\mathrm{I}}$, we can assign inversion-parity labels to the states at this point. Ab initio results indicate that both $|\varphi_{v_{2}\mathbf{k}}\rangle$ and $|\varphi_{c_{1}\mathbf{k}}\rangle$ are odd under $\hat{\mathrm{I}}$ exactly as $\mathbf{d}_{\Gamma\mathrm{K}}$, therefore the overall matrix element is zero.

We point out that the interband dipole is non-zero once the transition between a pair of bonding or antibonding combinations is considered. For example, taking into account the scattering $|\varphi_{v_{1}\mathbf{k}}\rangle\rightarrow|\varphi_{c_{1}\mathbf{k}}\rangle$, the intralayer-interband transition amplitudes sum constructively giving an intense overall interband dipole. The intense optical activity of these transitions is responsible for the main absorption peak denoted as C₁ in Fig. 3. Such prediction is confirmed by data in Fig. 6(b), where this quantity is shown as computed using the model (dashed blue line) and fully ab initio using Yambo (blue triangles). We notice that the reasonable agreement between the model and the ab initio results is an a posteriori confirmation of the validity of Eq. (9) to describe single particle states along $\Gamma\mathrm{M}^{\prime\prime}$.

As discussed in the Introduction, together with AB and AB^′ stackings, another stable bilayer C₃N motif is AA^′. Previous DFT calculations have effectively shown that these three stackings exhibit similar energies and coexistence of these motifs is expected at room temperatureWu et al. (2018). The crystal structure of AA^′-C₃N is shown in Fig. 7(a). This stacking has an inversion symmetry center (shown by the red dot in the figure), two mirror symmetry planes, represented by red dashed lines, and a two-fold rotation axis parallel to the stacking direction (indicated by the green dot). Interestingly, this stacking-motif is also invariant w.r.t. the non-symmorphic symmetry operation $\{\sigma_{xy}|\bm{\tau}\}$, corresponding to $z\rightarrow-z$ mirror symmetry followed by fractional translation of the vector $\bm{\tau}$, represented by the red arrow.

Structural optimization performed within PBE-D2 provides an in-plane lattice parameter $a$ = 4.849 Å, while the interlayer distance is equal to 3.22 Å, similarly to the other two stackings. For this relaxed atomic structure, the electronic structure at the DFT level is shown in Fig. 7(b). We notice that the highest valence state (found at the $\mathrm{M}^{\prime}$ point) has a higher energy than the lowest conduction state, occurring at the $\Gamma$ point, in agreement with the DFT results of Ref. [Wu et al., 2018] (our results give a negative ”gap” $\Gamma_{c}$ - M${}^{\prime}_{v}$ = -0.31 eV). Such metallicity is a problem related to the use of the PBE functional within the Kohn-Sham DFT scheme, since, experimentally, AA^′-C₃N has been shown to have a finite gap of about 0.4 eV Wenya et al. (2021) and hybrid DFT calculations Wenya et al. (2021) provide a semiconducting ground state.

The use of such metallic ground state to compute both QP corrections and optical properties is problematic, as it would induce a fictitious over-screening effect once the electron-electron screened interaction is evaluated using RPA approximation, providing inaccurate values for electronic gaps and exciton binding energies. Nevertheless, in the following, we will assume that the Kohn-Sham states computed at the PBE level are anyhow a good approximation for electronic wavefunctions, despite the problems of the associated Kohn-Sham energies.

Figure 7(b) indicates that the lowest direct band gap occurs almost in the middle of the $\Gamma\mathrm{M}^{\prime}$ direction. As a consequence, it is reasonable to expect that the low energy transitions which contribute to the lowest energy excitons also come from this portion of the BZ. Notably, $\mathbf{k}$ points along this direction are invariant both under the $\hat{\sigma}_{\Gamma\mathrm{M}^{\prime}}$ mirror symmetry and $\{\sigma_{xy}|\bm{\tau}\}$, therefore the electronic states can be properly labelled according to how they transform under these operations. Our DFT results indicate that the highest valence band is even under $\{\sigma_{xy}|\bm{\tau}\}$, while the lowest conduction is odd. In Fig. 7(b) these two states are indicated as $v^{+}_{2}$ and $c^{-}_{1}$, respectively. As the dipole operator $\mathbf{d}_{\epsilon}$ is invariant under $\{\sigma_{xy}|\bm{\tau}\}$ (assuming, as usual, incoming light with polarization direction $\hat{\epsilon}$ orthogonal to the stacking direction $z$), the matrix element $\langle\varphi_{v_{2}\mathbf{k}}|\hat{\epsilon}\cdot\mathbf{r}|\varphi_{c_{1}\mathbf{k}}\rangle$ is zero, independently of the direction of the polarization versor $\hat{\epsilon}$. Therefore, we expect light-induced scattering between these bands to be forbidden by symmetry and consequently the low-energy excitons composed by these transitions to be optically dark.

To confirm this symmetry-based analysis, we solve BSE computing the static electron-electron interaction in the direct kernel using PBE single particle wavefunctions and applying a rigid scissor $s^{W}_{0}$ to all the unoccupied bands, to manually remove DFT-spurious metallicity. Such scissor parameter has been chosen so that the minimum gap $\Gamma_{c}$ - M${}^{\prime}_{v}$ of the resulting band-structure was positive but smaller than the one found experimentally, to avoid under-screening effects. Furthermore, in the independent-particle part of the excitonic Hamiltonian, Eq. (1), we mimic quasi-particle corrections via a scissor operator applied to the DFT bands, manually chosen to obtain a minimum indirect band gap $\Gamma_{c}$ - M${}^{\prime}_{v}$ equal to the experimental one (0.4 eV). We underline that such scissor is therefore larger than the one introduced in the calculation of electronic screening.

The results of these calculations are shown in Fig. 8, where we have chosen $s^{W}_{0}=0.35$ eV , which fulfills the above-mentioned requirement to induce a small, positive band gap. The obtained spectra confirm the symmetry-based discussion just outlined. In fact, in both cases, low energy excitons (whose positions are indicated by black vertical bars in the inset) are optically dark independently of the polarization direction, and no absorption structure is observed in the low-energy region between 0.75 eV and 1.25 eV, where low-lying discrete excitons are found. The dark bound excitons are due to single particle transitions between bands $v^{+}_{2}$ and $c^{-}_{1}$ along the $\Gamma$-M^′ direction. This point is further clarified in the Supplemental Material sup , where, for completeness, we report the exciton wavefunctions of the low-lying excitons in reciprocal space. We point out that the observed optical quenching is stable with respect to small variations of the $s^{W}_{0}$ scissor parameter. In the Supplemental Material sup this is further confirmed by solving BSE using a RPA screening obtained with $s^{W}_{0}=$0.45 eV.

We emphasize that a more complete analysis of the optical properties of AA^′ stacking would require a better starting point for G₀W₀ and BSE calculations (i.e. a semiconducting electronic ground state). This is presently beyond the scope of this work and is left for future investigations.

In this work we discuss electronic and optical properties of bilayer C₃N using state-of-the-art ab initio calculations. As a first point, we find that optical absorption around 1.7 eV is favoured for all stacking structures studied. Furthermore, our results suggest a particular behaviour of BL-C₃N, i.e. the absence of low energy absorption peaks due to strongly bound excitons composed by electron-hole transitions with energies close to the minimum electronic direct gap.

These findings are explained in terms of independent particle effects as due to the negligible interband dipole between the lowest conduction and the topmost valence bands involved in the formation of the involved excitons. In the case of AB and AB^′ stackings, we develop a model for the single particle states of interest, to demonstrate that the overall interband dipole assumes negligible values because of the destructive interference of the contributions coming from the two layers. Furthermore, we also justify the quenching of low-energy absorption in AA^′-C₃N stacking, exploiting the symmetry properties of the crystal.

This work paves the way to future theoretical and experimental investigations on multilayer C₃N. On the one hand, it could be fascinating to investigate how optical properties of C₃N can be tuned varying the number of stacked monolayers or changing the twist angle between them. On the other hand, the abundant presence of dark or quasi-dark low-energy excitons in BL-C₃N could have important effects on exciton lifetimes and dynamics.

M.Z. thanks Marco Gibertini and Alberto Guandalini for stimulating discussions during this work and Nicola Spallanzani and Claudia Cardoso for support in the compilation of ab initio codes on different clusters. Lattice structures have been displayed using the XCrySDen Kokalj (1999) program. This work was partially funded by MaX – MAterials design at the eXascale – a European Centre of Excellence funded by the European Union’s program HORIZON-EUROHPC-JU-2021-COE-01 (Grant No. 101093374). It was also supported by the Italian national program PRIN2017 No. 2017BZPKSZ ”Excitonic insulator in two-dimensional long-range interacting systems” (EXC-INS).

We acknowledge ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing and ECOSISTER - Ecosystem for Sustainable Transition in Emilia-Romagna, funded under the National Recovery and Resilience Plan (NRRP) - NextGenerationEU - Mission 04, Component 2, Investments 1.4 and 1.5 [Call for tender n. 3277 dated 30/12/2021. Award Number:0001052 dated 23/06/2022].

M.J.C. acknowledges Brazilian agency CNPq and INEO (Nat. Inst. Sci. Technol. Organic Electronics - CNPq, FAPESP). Computational time on the Marconi100 machine at CINECA was provided by the Italian ISCRA program. We acknowledge EuroHPC Joint Undertaking for awarding us access to MeluXina at LuxProvide, Luxembourg. The authors gratefully acknowledge the LuxProvide teams for their expert support.