1 Introduction

While porphyrin and phthalocyanine have electronic structures that are markedly different [1,2,3,4,5], they appear as structurally analogous macrocyclic compounds, each comprising four pyrrole-like subunits interconnected to form a 16-member inner ring [6,7,8]. These macrocycles exhibit pronounced absorption bands in the visible and have remarkable thermal and chemical stability [9, 10]. Importantly, their photophysical and redox properties can be finely tuned by selecting appropriate metal centers or peripheral substituents [11,12,13,14,15]. They have emerged as important components in molecular materials exhibiting compelling electronic and/or magnetic properties [16,17,18]. Incorporating porphyrin or phthalocyanine moieties into polymer networks significantly enhances the versatility of the resulting materials [19,20,21,22,23,24]. Covalent organic frameworks (COFs) based on these macrocycles, hereafter referred to as Por-COFs or Pc-COFs, provide for a broad range of applications from catalysis, adsorption, and separation processes to therapeutic uses [25,26,27,28]. The properties leading to specific applications can be optimized via preparation strategies, structural designs, and monomer functionalizations [9,10,11, 15, 29,30,31,32,33,34].

In addition to the applications stemming from their ultrahigh porosity, Por- or Pc-COFs with a tetragonal lattice symmetry are of significant interest in the framework of organic electronics and spintronics [35,36,37,38,39,40,41,42]. Two-dimensional (2D) COFs with Por or Pc four-arm cores typically display a square lattice symmetry, such as a Lieb lattice or checkerboard lattice, as illustrated in Fig. 1. A Lieb lattice has a unit cell that contains one corner site and two edge sites and can lead to various exotic electronic properties derived from its intriguing combination of Dirac and flat bands [43,44,45]. The first realization of a Lieb lattice in COFs was reported in 2017 by Jiang and co-workers in the case of a pyrene-based COF [35]; chemical oxidation of this semiconducting material with iodine significantly enhanced its electrical conductivity, with the radicals appearing on the pyrene cores resulting in a high spin density and paramagnetism [35]. Subsequent computational studies have shown that these Lieb-like COFs can undergo magnetic phase transitions as the carrier doping concentration increases [46, 47]. In a recent theoretical work, Heine and co-workers [48] reported that several Zn-Pc COFs can have characteristic Lieb-lattice electronic bands, which can be shifted to the Fermi level via electron removal or atom substitution. The checkerboard lattice (which corresponds to the line graph of the square lattice) [49] has been studied as a lattice model for antiferromagnetism or flat-band superfluidity [50, 51], although its material design and synthesis have remained limited [36, 52]. The first demonstration of a checkerboard lattice was reported in 2018 by Crommie and co-workers in the case of a single layer of COF-420 [36]; COF-420 consists of porphyrin-based core units that are bridged by chemically asymmetric linkers. A topographic image obtained from scanning tunneling spectroscopy revealed its distinct checkerboard pattern [36].

Fig. 1
figure 1

Illustrations of typical lattice symmetries and representative cores and linkers for 2D tetragonal COFs. The four-arm cores and linkers are colored in green and blue, respectively. The grey square box in each lattice indicates the unit cell

In the case of hexagonal COFs, based on the determination of the frontier molecular-orbital (MO) symmetry of the building blocks (core and linker units) and the lattice symmetry, we described recently how one can tailor the electronic bands (Dirac band vs. flat band) near the Fermi level by selecting appropriate molecular units [53]. Here, we expand this investigation to tetragonal COFs, with the aim of providing a comprehensive understanding of their electronic structures. First, the generic band structures related to the square lattice, Lieb lattice, and checkerboard lattice are examined using tight-binding (TB) models, since these models offer a theoretical framework for analyzing the electronic properties of the various lattices. Then, through an analysis of the frontier MOs of the building units of several representative Por- or Pc-COFs and with the assistance of the TB models, we qualitatively predict the band structures of the corresponding COFs. In the following step, density functional theory (DFT) calculations are employed to validate these TB predictions and demonstrate the realization of the anticipated generic band structures in the 2D polymer networks.

2 Computational approaches

The ground-state geometry optimizations and electronic-structure calculations of the porphyrin, Zn-porphyrin (Zn-Por), phthalocyanine, and Zn-phthalocyanine (Zn-Pc) molecules were performed at the DFT level with the range-separated ωB97XD functional [54, 55] and the 6-31G** basis set (Gaussian16 package) [56]. The range-separation parameters (ω) were optimized nonempirically by minimizing \(J(\omega )={\left[{E}_{HOMO}(\omega )+IP(\omega )\right]}^{2}+{[{E}_{LUMO}(\omega )+EA(\omega )]}^{2}\), where \({E}_{HOMO[LUMO]}\) is the energy of the highest occupied molecular orbital (HOMO) [lowest unoccupied MO (LUMO)] and IP [EA] is the vertical ionization potential [electron affinity] [57]. The ω values were optimized to be 0.180, 0.178, 0.127, and 0.125 bohr-1 for the Por, Zn-Por, Pc, and Zn-Pc, respectively. Considering the frontier MOs of the molecules, we employed tight-binding models to investigate their generic band structures in the various tetragonal lattices. Further details on the tight-binding models are provided in the Supplementary Information (SI).

Turning to the molecular building units (cores and linkers) of the Por-/Pc-COFs, their ground-state geometries and frontier MOs were also evaluated at the DFT level using the ωB97XD functional [54, 55] and the 6-31G** basis set [56]. The range-separation parameter (ω) for each building unit was set to the same value as that tuned for the corresponding Por, Zn-Por, Pc, or Zn-Pc molecule. DFT band structure calculations on the 2D Por-/Pc-COFs were then performed at the level of the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) functional as implemented in the Vienna Ab initio Simulations Package (VASP) [58,59,60]. A Monkhorst-Pack k-points grid of 7 × 7 × 1 and a plane-wave cutoff of 500 eV were used. To ensure full decoupling between neighboring slabs, a vacuum layer over 15 Å was adopted. All atoms were allowed to relax until the atomic forces for structural relaxation were reduced to less than 0.01 eV/Å. The VASPKIT code [61] was used for post-processing of the VASP-calculated data and the VESTA program [62], to visualize the results.

3 Results and discussion

3.1 Frontier molecular orbitals of the porphyrin and phthalocyanine molecules

The frontier MOs of the porphyrin and phthalocyanine molecules (which exhibit D2h symmetry), are depicted in Fig. 2 along with those of their zinc-coordinated counterparts (D4h symmetry). The figure also collects the MO energies and symmetries. It is worth noting that: (i) The Zn-coordinated molecules have a doubly degenerate lowest unoccupied MO (LUMO), which in the case of the lower-symmetry metal-free molecules splits into the energetically close LUMO and LUMO + 1 levels. In all four instances, the next unoccupied MO is energetically well separated. Similarly, the highest occupied molecular orbital (HOMO) and the HOMO-1 levels are well separated with an energy separation over 1.5 eV. (ii) In Zn-Por and Zn-Pc, zinc does not contribute to the wavefunctions of the frontier levels that are entirely determined by the π-conjugated macrocycle.

Fig. 2
figure 2

Frontier molecular orbitals of the Por, Zn-Por, Pc, and Zn-Pc molecules, as calculated at the ωB97XD level of theory with the 6-31G** basis set. The MO energies and symmetries are also given. The black arrows indicate our convention for the x and y axes. (The wavefunctions of the HOMO-1 and LUMO + 2 levels are omitted as these levels are not considered in any of the subsequent calculations)

Building on the insight gained from previous studies [37, 41, 42, 53, 63, 64], these frontier MOs provide the orbital basis to initially elucidate the generic band dispersions in tetragonal Por- or Pc-COFs. In the next section, we illustrate various models where we consider one, two, or three orbitals per site. The relevance of these models to actual COFs will then be discussed later. We recall that, as was the case in our earlier work [63], the TB models look at the valence bands (VBs) and conduction bands (CBs) separately and can be different for the two cases.

3.2 Tight-binding model – square lattice

We begin with the simplest scenario, involving a single-orbital hopping (e.g., the single HOMO depicted in Fig. 2) within a square lattice, see Fig. 3a. Thus, here, there is only one site per unit cell, each site carrying one orbital. As a result, this lattice model yields a single band, as shown in Fig. 3b. The eigenvalue for the simple square lattice model is \(E\left(\overrightarrow{k}\right)=2t(\text{cos}{k}_{1}+ \text{cos}{k}_{2})\); \(k_{1,2}=\overrightarrow k\cdot\overrightarrow{a_{1,2}}\), and \(\overrightarrow{{a}_{1}}=\widehat{x}\), \(\overrightarrow{{a}_{2}}=\widehat{y}\). The bandwidth corresponds to 8|t|, i.e., eight times the electronic coupling strength between sites; this feature is commonly observed in polymer networks with tetragonal symmetry [37, 42, 65]. We note that the band is flat along the Y-X k-path (\(\left|{k}_{1}\pm {k}_{2}\right|=\pi\)), as illustrated in Fig. S1, which fulfills the phase cancellation condition (indeed, since \(\text{cos}{k}_{2}=-\text{cos}{k}_{1}\), this leads to \(E\left(\overrightarrow{k}\right)=0\)) [66].

Fig. 3
figure 3

Illustrations of (a) single-orbital hopping, (c) (πx, πy)-orbital hopping, and (e) (σ, πx, πy)-orbital hopping in a square lattice. The hopping integrals are indicated by arrows. The dashed square indicates the unit cell. The corresponding band structures obtained from TB models with (b) \({\varepsilon }_{0} = 0\) and \(t > 0\); (d) επ = 0, ppσ > 0, and ppπ = 0; and (f) εσ = -4.5|ssσ|, επ = 2.5|ssσ|, ssσ < 0, spσ =|ssσ|, ppσ =|ssσ|, and ppπ = 0. The inset in (b) depicts the first Brillouin zone. The ppσ, ppπ, ssσ and spσ hopping integrals are indicated by the arrows in (c) and (e)

For the next model, we consider the case where two frontier πx and πy molecular orbitals (MOs) need to be considered as the orbital basis (see Fig. 3c); these can be, for example, the degenerate LUMOs of Zn-Por or Zn-Pc. Here, two bands emerge with eigenvalues \({E}_{1[2]}=2*(\text{pp}\sigma *\text{cos}{k}_{1[2]}+\text{pp}\pi *\text{cos}{k}_{2[1]})\). Under the condition \(\left|{k}_{1}\right|=\left|{k}_{2}\right|\) along Γ-M1[2], the two bands are degenerate along the Γ-M k-path, as depicted in Fig. 3d. In the case of a negligible \(pp\pi\) hopping integral, the bandwidth corresponds to 4| \(pp\sigma\)|. We note that, when the core unit possesses only two-fold symmetry such as in Por or Pc and we need to consider the energetically close LUMO and LUMO + 1, since the degeneracy of these π orbitals is lifted, there appears an energetic splitting of the two bands along the Γ-M k-path (see Fig. S2c).

By incorporating an additional σ orbital, we can extend the orbital basis to (σ, πx, πy), as illustrated in Fig. 3e. This can be viewed as a combination of the earlier two lattice models with a single orbital and degenerate π orbitals and would have to be used, for instance, if a doubly degenerate LUMO is energetically close to the LUMO + 1. Then, the relative on-site energies of the σ and π orbitals, together with the differences in hopping integrals, can result in a variety of band dispersions, such as the one shown in Fig. 3f.

Given that the HOMO and LUMO levels of the Por and Pc systems exhibit opposite orbital parities (where orbital parity refers to the evolution of the wavefunction under spatial inversion, with even parity keeping the wavefunction the same upon inversion and odd parity changing the sign of the wavefunction), they provide the essential foundation for the emergence of topological states due to band inversion (see Fig. S3) [39, 41]. The (σ, πx, πy)-orbital basis in COFs is analogous to the (s, px, py)-orbital basis in inorganic materials, such as the Au/GaAs(111) interfacial system (even though the latter corresponds to a trigonal lattice) [67]. In Au/GaAs(111), the s- and p-orbital components of the three bands near the Fermi level were described effectively via an sp2 hybridization involving Au-s and As-p orbitals. Given the substantial spin-orbit coupling (SOC) effect in both Au and GaAs, the degeneracy at the Dirac point in the band structure is lifted, resulting in the formation of a topological insulating gap. Details regarding topological features in a square lattice with a (σ, πx, πy)-orbital basis are provided in the SI.

3.3 Tight-binding model – Lieb lattice

The Lieb lattice is a variant of the square lattice that consists of three sites within each square unit cell, as illustrated in Fig. 4a. Two of the sites (highlighted in blue) are connected to two neighboring sites, while the third site (marked in green) has four neighbors. Hereafter, these sites will be referred to as edge (blue) and corner (green) sites, respectively. When considering hopping only between nearest-neighbor (NN) sites, this geometric arrangement gives rise to an electronic band structure characterized by two distinctive features: (i) two dispersive bands forming a Dirac cone at the M point in the first Brillouin zone and (ii) a flat band intersecting the Dirac point (Fig. 4b). The flat band comprises electronic states localized exclusively on edge sites, whereas all sites contribute to the dispersive bands.

Fig. 4
figure 4

a Illustration of a single-orbital hopping in the Lieb lattice. The hopping integral is indicated by the black arrow. The spin-orbit coupling (SOC) term is considered as an imaginary hopping between the second nearest-neighbor sites, indicated by red arrows; \({\lambda }_{SOC}\) represents the magnitude of the SOC. The corner and edge sites are colored in green and blue, respectively. Examples of band structures obtained from the Lieb lattice models: (b) \({\varepsilon }_{corner/edge}\) = 0 and \({\lambda }_{SOC}\) = 0; (c) εcorner/edge = 0 and \({\lambda }_{SOC}\) = 0.1|t|; (d) \({\varepsilon }_{corner}\)= 0, \({\varepsilon }_{edge}\)= 0.2|t|, and \({\lambda }_{SOC}\) = 0; (e) \({\varepsilon }_{corner}\)= 0.2|t|, \({\varepsilon }_{edge}\)= 0, and \({\lambda }_{SOC}\) = 0; and (f) \({\varepsilon }_{corner}\)= 0.2|t|, \({\varepsilon }_{edge}\)= 0, and \({\lambda }_{SOC}\) = 0.1|t|

A SOC term can be incorporated into the Hamiltonian by considering an imaginary hopping between the second nearest-neighbor (2NN) sites (see Fig. 4a), in the spirit of the SOC term found in the Kane-Mele model [68]. The SOC effects can generally lift the degeneracy at the Dirac point, as illustrated in Fig. 4c and result in the opening of what are topologically non-trivial bandgaps above and below the flat band.

We recall that a topologically non-trivial bandgap is characterized by a topological invariant. In systems with inversion symmetry, the topological invariant known as the Z2 number [69] can be determined based on the parity at time-reversal invariant momenta [70]. A system with a nontrivial topology (Z2 = 1) represents a topological insulating phase that exhibits the quantum spin Hall effect. This phase is characterized by edge states with quantized conductivity. Additionally, the Z2 number is equivalent to the spin Chern number, \({C}_{s}=\frac{1}{2}({C}_{\uparrow }-{C}_{\downarrow })\). The Chern number is defined as: [71] \(C=\frac{1}{2\pi }\underset{BZ}{\int }{d}^{2}k{F}_{12}(k)\), where \({F}_{12}\left(k\right)=\frac{\partial }{\partial {k}_{1}}{A}_{2}\left(k\right)-\frac{\partial }{\partial {k}_{2}}{A}_{1}\left(k\right)\) is the Berry curvature; \({A}_{\mu }\left(k\right)=-i\langle {n}_{k}|\frac{\partial }{\partial {k}_{\mu }}|{n}_{k}\rangle\), the Berry connection; and \(|{n}_{k}\rangle\), the normalized wave function of the respective band.

We also examined scenarios with different on-site energies for edge and corner sites, which reflects the differences present between the COF core and linker units. In this context, the degeneracy of the Dirac and flat bands is disrupted, leading to a trivial bandgap characterized by a zero topological invariant, as depicted in Fig. 4d and e. The position of the bandgap, either below (\({\varepsilon }_{corner}\)< \({\varepsilon }_{edge}\)) or above (\({\varepsilon }_{corner}\)> \({\varepsilon }_{edge}\)) the flat band, is determined by the relative on-site energies between the two types of sites. Near the M point, the wavefunction components of these Dirac and flat bands are found to be markedly distinct from those of the ideal Lieb lattice model (see Fig. 4b). Specifically, the wavefunction at the M point of the separated Dirac band (the bottommost band in Fig. 4d or the topmost band in Fig. 4e) primarily comes from the corner sites, while the edge sites predominantly contribute to the flat band and the other Dirac band.

We also considered a case where a SOC term is incorporated in addition to the on-site energy difference. The resulting band structure has three isolated bands, as depicted in Fig. 4f; the two bandgaps at the M point correspond to a topological insulating phase with a nonzero Z2 number. By adjusting the parameters of the on-site energy difference and the SOC strength, the Chern number of the middle flat band can be either zero or nonzero [46].

Given the four-fold symmetry of the corner site in a Lieb lattice, the doubly degenerate frontier MOs of the four-arm core units can also emerge in tetragonal COFs. We examined a special case of the Lieb lattice where the corner site carries degenerate (πx, πy)-orbitals and the edge site contains a single orbital, as depicted in Fig. 5a. By considering only the NN interactions between corner and edge sites, while excluding interactions between corner-to-corner or edge-to-edge sites, the band structure obtained from the TB model is presented in Fig. 5b. This band structure is markedly different from those shown in Fig. 4. Therefore, the relative on-site energies for the corner and edge sites as well as the MO degeneracy of the corner site play an important role in determining the electronic band structure of the Lieb lattice.

Fig. 5
figure 5

a Illustration of (σ, πx, πy)-orbital hopping in the Lieb lattice. The corner site carries degenerate (πx, πy)-orbitals and the edge site, a single orbital. Only nearest-neighbor interactions are considered; the hopping integral spσ is indicated by the arrows between the corner and edge sites. b Band structures obtained from this lattice model: εσ = 0.5|spσ|, επ = -0.5|spσ|, and spσ > 0

3.4 Tight-binding model – checkerboard lattice

A checkerboard lattice is created by introducing additional sites at the center of each square in the basic square lattice, as depicted in Fig. 6a. When considering only NN hopping, the TB model in the case of a single orbital per site yields the two bands shown in Fig. 6b; the bands have a folded structure derived from those shown in Fig. 3b. If the 2NN hopping parameter is non-zero, the partial flat band at half filling gains some dispersion, as illustrated in Fig. S4a. By extending the orbital basis to include two (πx, πy) or three (σ, πx, πy) orbitals per site, the typical band structures that can be obtained are depicted in Fig. 6c and d, respectively. Similarly, the partial flat bands along the M-Γ k-path gain some dispersion when considering a nonzero ppπ parameter (see Fig. S4b).

Fig. 6
figure 6

a Illustration of the single-orbital hopping model in the checkerboard lattice; the nearest-neighbor hopping integral is indicated by t1 and the second nearest-neighbor hopping, by t2. Band structures obtained from the checkerboard lattice models: (b) single-orbital hopping with nearest-neighbor hopping integral \(t_{1} > 0\) and \({\varepsilon }_{0} = 0\); (c) (πx, πy)-orbital hopping with επ = 0, ppσ > 0, and ppπ = 0; (d) (σ, πx, πy)-orbital hopping with εσ = -4.5|ssσ|, επ = 2.5|ssσ|, ssσ < 0, spσ =|ssσ|, ppσ =|ssσ|, and ppπ = 0

It is worth noting that the band structure derived from the checkerboard lattice model exhibits similarities with those of materials possessing a Cairo pentagonal lattice symmetry. In each 2D pentagonal lattice (see Fig. S5), there are two types of sites: one is fourfold-coordinated and the other is threefold-coordinated, with a total of six sites per unit cell. The fourfold-coordinated sites display the checkerboard lattice symmetry. Considering one orbital per site, this lattice model yields six bands. For reference, the checkerboard lattice with a (σ, πx, πy)-orbital basis also produces six bands, which indeed look like those obtained from the pentagonal lattice model [40]. Since the theoretical prediction of penta-graphene in 2015 by Zhang and co-workers [72], a variety of other 2D pentagonal materials have been theoretically predicted and/or experimentally synthesized [73,74,75,76,77]. These materials are reported to exhibit tunable properties ranging from semiconducting to metallic and even topologically insulating states [78].

3.5 Covalent organic frameworks

We selected three COFs with either Zn-Por or Zn-Pc as the four-arm core units to serve as representative examples. These examples were chosen to demonstrate that the characteristic band dispersions associated with the square lattice, the Lieb lattice, and the checkerboard lattice, which we described above, can be realized in 2D polymer networks.

3.5.1 Square lattice

The first COF we examined is built from zinc-5,10,15,20-tetraethynylporphyrin (Zn-TEP) units (we refer to it as Zn-TEP COF) [38]. Its chemical structure is illustrated in Fig. 7a. Prior to conducting band-structure calculations, we first evaluated the frontier MOs of the Zn-TEP unit, as shown in Fig. S6. Based on the understanding gained from the TB models we discussed above, we can predict the characteristics of the bands appearing near the Fermi level. Following the square lattice model with single-orbital hopping (see Fig. 3b), the topmost VB of the Zn-TEP COF is anticipated to be a single band derived from the Zn-TEP HOMO. The two lowest CBs are expected to originate from the degenerate LUMO levels of Zn-TEP, following the (πx, πy)-orbital hopping model in a square lattice, which results in the typical band dispersions illustrated in Fig. 3d.

Fig. 7
figure 7

Two-dimensional Zn-TEP COF: a Chemical structure; the red dashed box represents the unit cell. b DFT-PBE-calculated band structure; the VBs and CBs are colored in blue and red, respectively. c Illustration of the wavefunctions at the Γ point for bands 1 to 4, as labelled in (b)

The DFT-calculated band structure of Zn-TEP is presented in Fig. 7b and confirms our predictions: The lowest two CBs display the characteristic band dispersion arising from a square lattice model with degenerate (πx, πy)-orbital hopping, while the top VB is solely derived from the HOMO of the Zn-TEP. Figure 7c displays the wavefunctions of bands 1 to 4 at the Γ point. We note that band 3 is a flat band, originating from the HOMO-1 level of Zn-TEP; this flatness is related to the fact that the acetylene units do not contribute to the HOMO-1 wavefunction (see Fig. S6) [42]. A comparison of the Γ-point wavefunctions with the frontier MOs of Zn-TEP (Fig. 7c vs. Fig. S6) corroborates that the electronic bands near the Fermi level are indeed derived from the frontier MOs of the corresponding building unit. This confirms that, by examining the characteristics of the frontier MOs of the building units along with the COF lattice symmetry, one can effectively predict the band characteristics in the corresponding COF with the assistance of TB models.

3.5.2 Lieb lattice

We recall that the Lieb lattice has two types of sites: corner sites and edge sites. For the interaction between these sites to be effective, it is essential to have a close alignment of their on-site energies; this facilitates the desired electronic coupling between the two types of sites and leads to the emergence of electronic bands characteristic of the Lieb lattice. Therefore, to realize the electronic features of a Lieb lattice in a 2D COF, it is important to select a linker unit with an on-site energy comparable to that of the core unit. The Zn-Pc-based COF with pyrazine-pyrene-pyrazine linkers [79] illustrated in Fig. 8a, is a pertinent example in this regard. The frontier MOs of the building units are shown in Fig. S7. Given that the Zn-Pc HOMO level is much higher in energy than that of the pyrene-pyrazine linker (by over 1 eV), the top VB is expected to be exclusively derived from the former level. Subsequently, the appropriate model for the top VB in fact corresponds to a single orbital hopping in a square lattice, with the anticipation that the band will be flat as only the corner sites will contribute to it. Turning to the unoccupied levels, the degenerate Zn-Pc LUMOs levels are lower in energy than the linker LUMO but with a smaller energetic separation of approximately 0.7 eV; thus, the bottom CBs of this COF can be expected to display a band structure akin to that illustrated in Fig. 5b. In addition, since the LUMO + 1 level of Zn-Pc and the LUMO + 2 level of the linker are rather closely aligned (with a separation of about 0.48 eV), the electronic states derived from these orbitals are anticipated to exhibit Lieb-like characteristics.

Fig. 8
figure 8

Two-dimensional Zn-Pc-based pyrazine-linked COF: a Chemical structure; the red dashed box represents the unit cell. b DFT-PBE-calculated band structure; the VBs and CBs are colored in blue and red, respectively. c Illustration of the wavefunctions at the Γ point for bands 1 to 3, as labelled in (b)

The DFT-calculated band structure is shown in Fig. 8b. Indeed, the dispersions of the upper unoccupied bands in the energy window between 2 and 3.5 eV for the CBs exhibit three bands (denoted 1 to 3) possessing the characteristics of a Lieb lattice. An analysis of their wavefunctions at the Γ point, shown in Fig. 8c, reveals that bands 1 and 3 originate from both the core and linker units, while band 2 is predominantly contributed by the LUMO + 2 level of the linker units; these attributes are consistent with those depicted in Fig. 4e.

The wavefunctions of the bands near the Fermi level are depicted in Fig. S8. As predicted, the topmost VB consists of a single band originating solely from the Zn-Pc HOMO. The bottom CBs can be classified into three groups (CB, CB + 1, and CB + 2), each containing two bands. There exists a close correspondence between these CBs and the Zn-Pc LUMO and the LUMO and LUMO + 1 levels of the pyrazine-pyrene-pyrazine linker (see Fig. S7). In contrast to the width of approximately 1 eV seen for the Lieb bands, the widths of the CBs and VBs near the Fermi level are much narrower and smaller than 0.1 eV. This can be attributed to the very limited electronic coupling between the pyrazine-pyrene-pyrazine-based linkers and Zn-Pc units. By modulating the degree of conjugation in the linker unit (i.e., by varying the linker length), the relative positions of the three bands derived from the Lieb lattice can be adjusted to be in closer proximity to the Fermi level [48]. Again, in the framework of the COF lattice symmetry, the analysis of the symmetries and energies of the frontier MOs of the core and linker units allows for the prediction of the band dispersion near the Fermi level in a specific COF.

3.5.3 Checkerboard lattice

When the core units situated at both the center and corner sites are identical, the primitive unit cell of the COF adopts the structure of a simple square lattice. However, introducing a relative orientation between these two sites, as exemplified by the Zn-Pc-based anthracene-linked COF [40] depicted in Fig. 9a, leads to a unit cell with the symmetry of a checkerboard lattice. The frontier MOs of the Zn-Pc core and the tribenzo[f,k,m]tetraphene linker are shown in Fig. S9. The degenerate Zn-Pc LUMO levels are significantly lower in energy than the linker LUMO (difference of ca. 1.6 eV), while the Zn-Pc HOMO is higher than that of the linker by about 1 eV. Thus, the bottom CBs and top VBs of this ZnPc-based anthracene-linked COF are expected to display band dispersions similar to those depicted in Fig. 6.

Fig. 9
figure 9

Two-dimensional ZnPc-based anthracene-linked COF: a Chemical structure; the red dashed box represents the unit cell. b DFT-PBE-calculated band structure; the VBs and CBs are colored in blue and red, respectively. c Illustration of the wavefunctions at the Γ point for bands 1 to 6, as labelled in (b)

The band structure calculated at the DFT-PBE level of theory is shown in Fig. 9b and indeed exhibits features like those in Fig. 6. Figure 9c presents the wavefunctions corresponding to bands 1 to 6 at the Γ point and confirms that the bottom CBs originate from the degenerate LUMOs of Zn-Pc, while the top VBs are derived from the single HOMO of Zn-Pc. The wavefunctions at the Γ point for bands 7 and 8, depicted in Fig. S10, clearly carry the orbital characteristics of the HOMO level of the tribenzo[f,k,m]tetraphene linker. These results underline again that the most relevant characteristics of the band structure of a specific COF can be qualitatively predicted by analyzing the symmetry and energy of the frontier MOs of its building units in conjunction with the lattice symmetry.

In the case of Por-/Pc-based COFs with two-fold symmetric porphyrin or phthalocyanine as four-arm core units, the band structures closely resemble those of COFs with Zn-Por/Pc as core units [41, 80], despite the lifted degeneracy at the M point due to the energy splitting between the LUMO and LUMO + 1 levels in the metal-free porphyrin and phthalocyanine molecules (Fig. 2). The comparison between Zn-Pc-based and Pc-based COFs is presented in Figs. S11 and S12.

While our present work focuses on the electronic structure of COF monolayers, we have also briefly investigated the impact of the interlayer stacking patterns on the electronic properties by using Zn-TEP COF as a representative example, as detailed in the SI. Due to the relatively small bandgap of these Pc-based COFs, they are promising materials to achieve band inversion and a topological state. A topological phase transition from a conventional insulating state to a topological Dirac semimetal state has been theoretically predicted to occur in a Pc-based COF upon the application of external strain [41]. However, it is worth bearing in mind that the realization of a first-order topological insulating phase usually requires a significant SOC strength to facilitate the band inversion and bandgap opening. In this regard, it is thus challenging to achieve a topological insulating state in COFs composed solely of light elements (such as H, C, N, O, and F) or where only light elements contribute to the frontier wavefunctions. We note that incorporating heavy metal atoms with a strong SOC effect is an effective approach to realizing a topological insulating phase in porphyrin/phthalocyanine-based metal-organic frameworks (MOFs). For instance, phthalocyanine-based MOFs with metal centers other than zinc have been studied and shown to exhibit a topological insulating phase, leading to the quantum spin/anomalous Hall effect [39, 81].

4 Conclusions

Employing tight-binding models and density functional theory calculations, we have investigated the electronic structures of tetragonal COFs, using porphyrin- and phthalocyanine-based COFs as representative examples. The generic band structures corresponding to the square lattice, the Lieb lattice, and the checkerboard lattice were initially explored through tight-binding model analyses. The realization of these band structures was demonstrated in 2D polymer networks that incorporate porphyrin or phthalocyanine as building units, thereby confirming the versatility and applicability of these TB models in the context of COFs.

Expanding on our earlier investigations of hexagonal COFs, we have determined that the band structures in tetragonal COFs can also be tailored by choosing the symmetries of the core units and/or by adjusting the relative energies of the frontier MOs of the core and linker units. The approaches taken to engineer the band structures outlined in this and our previous studies are expected to be applicable to a wide variety of COFs, thereby broadening the ability to tune the electronic properties and enhancing the functionality of these materials.

Finally, we emphasize that, unlike the first-order topological insulating phase that requires a substantial spin-orbit coupling effect to achieve band inversion and bandgap opening, the topological semimetal phase and higher-order topological insulating phase can be realized in COFs composed exclusively of light elements [41, 63, 82].This distinction highlights the versatility of COFs in accommodating various topological phases, thereby expanding their potential applications in the realm of topological materials science.