Vicinal metal surfaces as potential catalysts for phosphorene epitaxial growth

TL;DR

This study computationally explores how stepped vicinal metal surfaces, like Cu(311), can favor the growth of black phosphorene (BLK-P) over blue phosphorene (BL-P), offering new avenues for phosphorene synthesis.

Contribution

It demonstrates that stepped metal surfaces can preferentially stabilize BLK-P, providing a novel approach for epitaxial growth of phosphorene.

Findings

Surface steps favor BLK-P over BL-P on Cu(311).

Surface density of states and orbital hybridization explain the stabilization.

Vicinal metal surfaces may enable phosphorene epitaxial growth.

Abstract

Phosphorene, a single layer of black phosphorous (BLK-P), has a significant potential for flexible and tunable electronics, but attempts to grow it epitaxially have been unsuccessful to date. Meanwhile, hexagonal blue phoshorous (BL-P) has been achieved on closed-packed (111) metal surfaces in special growth conditions of high vapor pressure and high reactivity of phosphorous. The (111) surfaces favors BL-P over BLK-P due to its hexagonal symmetry. Here, we investigate computationally the alternative offered by stepped substrates. Using the Cu(311) surface as a model, we find that surface steps can favor energetically BLK-P over BL-P. This can be rationalized in terms of surface density of states and orbital hybridization, which lead to a stronger surface bonding of the lower BLK-P half-layer. This work suggests that vicinal metal surfaces of metals can offer a viable path towards phosphorene synthesis.

Tables

Table 1. Table 1: Lattice mismatches δ 𝛿 \delta and calculated formation energy E f subscript 𝐸 𝑓 E_{f} (Eq. 1 ) and adsorption energy E a ​ d ​ s subscript 𝐸 𝑎 𝑑 𝑠 E_{ads} (Eq. 2 ) of BL-P and BLK-P on Cu ( 111 ) 111 (111) and Cu ( 311 ) 311 (311) substrates. The subscripts TS and D2 indicate the vdW correction scheme. δ 1 subscript 𝛿 1 \delta_{1} and δ 2 subscript 𝛿 2 \delta_{2} measure lattice mismatches along y 𝑦 y and x 𝑥 x directions, respectively. A negative value of δ 𝛿 \delta indicates that the overlayer is compressively stressed. All energy values are given in meV.

		Cu$\left(111\right)$	Cu$\left(311\right)$
BLK-P	$E_f$	-285	-386
	$E_{f-D2}$	-521	-653
	$E_{f-TS}$	-548	-695
	$E_{ads}$	-170	-271
	$E_{ads-D2}$	-316	-448
	$E_{ads-TS}$	-342	-489
	${\delta }_1$ (%)	-2.8	-7.1
	${\delta }_2$ (%)	-2.4	-2.6
BL-P	$E_f$	-348	-338
	$E_{f-D2}$	-632	-620
	$E_{f-TS}$	-667	-658
	$E_{ads}$	-235	-225
	$E_{ads-D2}$	-449	-437
	$E_{ads-TS}$	-485	-476
	${\delta }_1$ (%)	-1.6	+0.3
	${\delta }_2$ (%)	-1.3	-1.7

Equations

$$E_f=(E_{slab+P}-E_{slab}-N_P\times E_{P_4}/4)/N_P.$$

$$E_{ads}=(E_{slab+P}-E_{slab}-E_P)/N_P,$$

Full text

Mono or few-layer two-dimensional (2D) materials are extremely flexible, have atomically abrupt interfaces, and in principle can be assembled into 3D structures without typical lattice matching requirements characteristic to covalent epitaxy. The few-layer materials are held together by weak van der Waals (vdW) forces. In recent years BLK-P has emerged theoretically as a promising quasi (2D) electronic material characterised by high mobility, ambipolar transport and tunable bandgap Liu et al. (2014). Bulk black phosphorous is available for experiments, which inspired exfoliation of BLK-P by analogy to graphene exfoliated from bulk graphite. Exfoliated multiple-layer BLK-P flakes has remained the only viable method to obtain BLK-P. In contrast to other 2D materials which were grown on metal surfaces (e.g., graphene Loginova et al. (2009), hexagonal boron nitride (hBN) Siegel et al. (2017), silicene Vogt et al. (2012), germanene Cantero et al. (2017)), the epitaxial growth of BLK-P was not succesful to date, limiting its lateral size and alignment with various substrates. While the growth of 2D materials mentioned above on $(111)$ surfaces was enabled by their $C_{3v}$ symmetry, the $C_{2v}$ symmetry and strong in-plane anisotropy of BLK-P preclude the growth of BLK-P on these substrates. Instead, the theoretically predicted hexagonal BL-P Zhu and Tománek (2014a) has been obtained Zhang et al. (2016).

We attempted to develop a thin film growth epitaxiual method for BLK-P by using the chemical vapor deposition (CVD) applied previously with success to the growth of (layered) hBN Siegel et al. (2017). We used either phosphane PH3 or tributyl phosphate TBP precursors. Even at temperatures as low as $300$°C on copper surfaces we were able to obtain only metal phosphides (P dissolved in the substrate). This is inferred from X-ray diffraction and Raman characterisation, which show a very clear structural change from polycrystalline Cu to Cu3P. Similarly, we have observed formation of phosphides on Ni, Au, and Ag at low temperatures. This highlights the high reactivity of P precursor fragments under CVD conditions. It confirms the conclusion of modeling efforts suggesting that flakes of BLK-P are not structurally stable on low-index transition metal surfaces Gao et al. (2016). It may be possible that the growth of BLK-P be achieved by molecular beam epitaxy (MBE), as in the case of BL-P. For that, an appropriate metal substrate needs to be found that would favor the formation of BLK-P over BL-P. Several routes can be followed to alter the catalytic activity of metals, including the demonstration of hybrid metal-hBN substrate with an ‘adjusted van-der Waals bonding’ that could stabilize the BLK-P Gao et al. (2016), tailoring the chemical activity of the substrate by alloying, or using stepped surfaces Davidson and Pickett (1994) with lower symmetries and increased chemical activity in one direction. In this study, we follow the latter route.

In this paper, we investigate computationally the possibility raised that stepped surfaces catalyse the growth of BLK-P. This approach is in the spirit of studies of BLK-P vs BL-P energetics on low index metals or semiconductors, which indicated that low-index surfaces favor BL-P over BLK-P, with the exception of a theoretical Sn surface Zeng et al. (2017); Qiu et al. (2018). To study the increased catalytic strength induced by surface step we used Cu since its low-index surface has a lower catalytic activity by comparison to, e.g., Ni, Ir or Fe Bokdam et al. (2014). We chose the Cu$(311)$ terrace as a minimal system size for periodic cell calculations that has also a relatively small lattice mismatch with BLK-P compared to other metals.

The calculations were performed using the Vienna ab initio simulation package (VASP) Kresse and Hafner (1993); [Kresse and Furthm

](#bib.bib14) and based on frozen-core full-potential projector augmented-wave method (PAW) [Bl

](#bib.bib15) with the Perdew, Burke and Ernzernhof Perdew et al. (1996) generalized gradient approximation. The convergence with respect to the size of $k$-sampling, cutoff radius, substrate and vacuum thickness was checked. These resulted in a $k$-point mesh of 8$\times$8$\times$1 for density of states (DOS), a plane-wave cutoff radius of $400$ eV, a substrate thickness of six Cu layers, and a vacuum distance between the periodic images of $20$ Å. Dispersion corrections to standard density functional theory must be included in modeling surface chemistry of 2D materials. Among the vdW -corrected functionals, Grimme’s empirical long range dispersive force correction (D2) Grimme (2006) has been able to qualitatively capture the charge redistribution and interlayer spacing of BLK-P system Shulenburger et al. (2015) and has been shown to be accurate for predicting the energetics of rare gas atoms on metal surfaces Chen et al. (2011, 2012). In this work, we used this correction during structural relaxation. We used the conjugate-gradient algorithm for structural relaxation with a force convergence criterion of $0.01$ eV/Å. For the relaxed structures, we used also the self-consistent method of Tkatchenkoand Scheffler (TS) Tkatchenko and Scheffler (2009) a posteriori to verify the ordering of converged energies.

The calculated lattice constants of freestanding BLK-P are a_{1(0)}^{BLK-P}$$=$$4.59 Å and a_{2(0)}^{BLK-P}$$=$$3.30 Å (armchair and zigzag directions, respectively). For freestanding BL-P, the lattice constant is a^{BL-P}$$=$$3.27 Å; for strain reference we use the zigzag and armchair lattice vector parameters a_{1(0)}^{BL-P}$$=$$5.66 Å and a_{2(0)}^{BL-P}$$=$$3.27 Å  (armchair and zigzag directions, Fig.1). The bulk metal lattice constant was estimated to be a^{Cu}$$=$$3.63 Å. The corresponding lattice mismatches along the principal directions of BLK-P or BL-P are listed in Tab. 1. The standard dipole correction for surface adsorption calculations was used to guarantee convergence. The clean metal slab was first relaxed with the three bottom layers fixed at their ideal bulk positions. The top three metal layers were further relaxed again in the presence of BLK-P or BL-P. For structural relaxations of the adsorbed system in the large supercells, we used only the $\Gamma$ point sampling of the Brillouin zone.

To match BLK-P on a finite Cu slab with a minimal strain [Fig.1(a)], the supercell contains four BLK-P unit cells with the zigzag direction aligned with five Cu$(311)$ surface unit cells along the step ($x$ direction), and four BLK-P unit cells with the armchair direction aligned with four Cu$(311)$ surface unit cells perpendicular to the step ($y$ direction). This corresponds to $240$ Cu atoms in the Cu$(311)$ slab and to N_{P}$$=$$64 P atoms in the BLK-P included in the supercell, or a 4:5 ratio of P:Cu atoms at the interface (in the bottom BLK-P layer and the topmost metal layer). The interface structure of BLK-P/Cu$(311)$ shows that all Cu atoms at the Cu$(311)$ surface steps are shared between the two rows of BLK-P atoms and form a bond, while the Cu atom on the terrace do not form bonds to BLK-P atoms. This is because the lattice constant $a_{1(0)}^{BLK-P}$ of freestanding BLK-P is $7\%$ larger than the distance of $4.26$ Å  between steps on Cu$(311)$. The lattice mismatches are listed in Table 1 and defined as $\delta_{i}$= $(a_{i}-a_{i(0)})/a_{i(0)}$, where ${i}$ indexes the Cartesian components. Center-to center Cu-P bond length ranges from $2.25$ Å to $2.43$ Å, with an average Cu-P bond length of $2.32$ Å. However, the averaged interfacial vertical distance ($d$) between the bottom of BLK-P and the top most Cu layer is $1.68$ Å. The average bond length above is close to the sum of covalent radii of Cu ($1.27$ Å) and P ($1.08$ Å) Zhu and Tománek (2014b), i.e., $2.35$ Å. This indicates the possibility of a considerable orbital overlap between P and Cu orbitals leading to chemical bonding.

In comparison to BLK-P, BL-P exhibits a different adsorption pattern on Cu$(311)$ [Fig.1(b)]. In this case, even though all the Cu step atoms form bonds with a single row of BL-P atoms, every other Cu atom on the terrace also from bonds to some of the BL-P atoms. The supercell contains four BL-P unit cells with the zigzag direction aligned with five Cu$(311)$ surface unit cells along the step, and three BL-P unit cells along the zigzag direction aligned with four Cu$(311)$ surface unit cells perpendicular to the step. The supercell includes $240$ Cu atoms and N_{P}$$=$$48 P atoms, or a 3:5 ratio of P:Cu atoms at the BL-P/metal interface. The average lattice constants of BL-P are slightly increased to a_{1}^{BL-P}$$=$$5.68 Å ($0.3\%$ increase) to match the Cu unit cell along the $y$ direction, and decreased to a_{2}^{BL-P}$$=$$3.22 Å ($1.7\%$ reduction) along the $x$ direction. Center-to center Cu-P bond length ranges from $2.26$ Å  to $2.50$ Å, with an average Cu-P bond length of $2.33$ Å. The optimized interfacial distance ($d$) between the bottom of BL-P and the top most Cu layer is found to be $1.90$ Å, which is $0.22$ Å  greater than BLK-P on Cu$(311)$.

A comparison between phosphorous allotropes bonding on metal substrate is based on the formation energy $E_{f}$. This is defined with respect to a common reference state for the energies of P atoms, chosen here to be the tetraphosphorus P4 molecule:

[TABLE]

Here, $E_{slab+P}$, $E_{slab}$, and $E_{P_{4}}$ are the total energies of the slab with phosphorene on it, of the metal slab, and of P4, respectively. This picture assumes that the allotropes are assembled from fragments of the P4 precursor dissociated on the surface. Another comparison can be made in terms of the adsorption energy $E_{ads}$, often used in the literature:

[TABLE]

where $E_{P}$ is the energy of the isolated phosphorene allotrope. This assumes an allotrope already formed before it is adsorbed on the surface, such as in an exfoliation/transfer process. The calculated $E_{f}$ and $E_{ads}$ are summarized in Tab. 1. Negative values indicate bonding. Both BL-P and BLK-P bind to either surface, with BL-P bonding more strongly than BLK-P on Cu$(111)$, and BLK-P bonding more strongly that BL-P on Cu$(311)$. We note that the difference between $E_{ads}^{BLK-P}$ and $E_{ads}^{BL-P}$ follows closely the difference E_{f}^{BLK-P}$$-$$E_{f}^{BL-P} for all chemistry levels considered. In addition, for both allotropes there is significant bonding already at the PBE level. The addition of dispersive force corrections further increase the bonding energies, maintaining the ordering BL-P, BLK-P. For BLK-P $E_{ads}$ on the Cu$(311)$ surface is relatively large, i.e., $0.49$ eV. Using the $E_{f}$ criterion, we obtain that BLK-P is favored by \approx$$40 meV/atom on Cu$(311)$, whereas using $E_{ads}$ the difference is reduced to \approx$$10 meV/atom. The former criterion seems more realistic from the point of view of assembling the P allotropes from P fragments on the surface and indicates that BLK-P is stable on Cu$(311)$ against transition to BL-P at room temperature.

An insight into the relation between geometries and adsorption energies can be obtained from the change in local electronic density at the interface defined as \Delta\rho(\bm{r})$$=$$\rho_{slab+P}$$-$$\rho_{slab}$$-$$\rho_{P}. Figs. 2(a),(b) show $\Delta\rho(\bm{r})$ for both allotropes using the same iso-values. The blue isosurfaces represent local charge migrating to new, hybrid orbitals formed by the P-Cu interaction. The red isosurfaces represent migration of electronic charge from the linear superposition of the separate Cu slab and P allotrope densities. In addition, Figs. 2(c,d) show the plane average \Delta\sigma_{\pm}(z)$$=$$\int\Delta\rho_{\pm}(\bm{r})dxdy for charge migration to hybrid orbitals ($+$, blue) and depletion of the linearly-superposed orbitals ($-$, red). The topmost Cu plane and the lowest P plane are represented by dotted lines. We notice that $\Delta\sigma_{-}$ peaks on these planes and are separated by a charge minimum. On the other hand, $\Delta\sigma_{+}(z)$ is maximum at the interface between these planes and decays towards the two components of the system, with satellite minima on the second layer of Cu atoms and the topmost layer of P atoms. This is consistent with the $\Delta\rho$ iso-surface plots. We interpret this as an overall bonding between P and Cu atoms via hybrid states. A larger population in the hybrid orbitals is associated with stronger bonds, reflected in the adsorption energy. In the case of BLK-P, $\int\Delta\sigma_{+}(z)dz$ gives $0.166~{}e$ per surface Cu atom of the Cu$(311)$ substrate. For BL-P, $\int\Delta\sigma_{+}(z)dz$ amounts to $0.146~{}e$ per Cu atom per surface Cu atom in Cu$(311)$. The net difference of $13.75$% in charge transfer per surface Cu atom favors BLK-P. On the other hand, in the case of the Cu(111) $\int\Delta\sigma_{+}(z)dz$ amounts to $0.135~{}e$ per surface Cu atom for BLK-P, and to $0.154~{}e$ per surface Cu atom for BL-P. The net difference is $12.60$% in favor of the latter. Besides the increased symmetry, a factor for the change in preference comes from the P:Cu ratios at the interface, i.e., 4:5 and 16:25 for BLK-P and BL-P on Cu(111). This points to a slightly increased number of interface P-Cu bonds for BL-P on Cu(111) in comparison to Cu(311). These results are consistent with the ordering of adsorption energies.

The discussion above on charge redistribution into hybrid states is a global view over all orbital types and energies. The nature of the hybrid states can be understood by analyzing either the partial density of electronic states (PDOS) or by a Crystal Orbital Hamilton Populations (COHP) analysis Deringer et al. (2011); Dronskowski and Bloechl (1993). The latter emerged from the Frontier Orbital Theory Hoffmann (1988) and is directly related to the PDOS by multiplication with bond orbital populations extracted from the first-principles calculation. Two features usually occur upon absorption: a hybridization of metal Fermi level and the adsorbate conduction band minimum (or lowest unoccupied orbital), which give rise to an occupied state below the Fermi level and to an empty orbital higher in the metal conduction band; and a strong splitting of the deep adsorbate states. The latter usually results from a hybridization of the upper part of the adsorbate valence band (or highest occupied orbital) with the first unoccupied metal state. In our case, PDOS show that BLK-P undergoes a larger splitting of the deep P orbitals than BL-P. The magnitude of the splitting is associated with a larger bonding energy and is consistent with symmetry arguments. Namely, the laterally-extended $p$ orbitals of BLK-P have a $C_{2v}$ coordination Hoffmann (1988) that prefers the symmetry of the $(311)$ surface, whereas the vertically-extended $p$ orbitals of BL-P have a $C_{3v}$ symmetry that prefer binding with vertical metal orbitals on Cu$(111)$. This is reflected in the COHP analysis which decomposes the bonding in shows the contributions from each orbital pair.

Figs. 3,4 show the COHP analyses on the $(311)$ surface. These are normalized by the number of bonds and are represented with a minus prefactor so that positive values of –COHP represent bonding, and negative values antibonding contributions. We show the contributions from the Cu${}_{d_{z^{2}}}$ and Cus only with all P orbitals, since we find that these have the largest –COHP peaks. In Fig. 3, we notice strong bonding contributions from Cu${}_{d_{z^{2}}}$ in the interval [$-6$,$-4$] eV for both allotropes, but also large antibonding contributions in the interval [$-3$,[math]] eV which cancel these bonding contributions. On the other hand, Fig. 4 shows strong bonding contributions from Cus-Ps in the interval [$-16$,$-8$] eV for both allotropes, accompanied by much smaller Cus-Ps antibonding contributions in the interval from [$-8$,[math]] eV. In addition, the Cus-Pp bonding contributions between [$-7$,$-4$] eV for BLK-P, and [$-7$,$-6$] eV for BL-P have very little cancellation from antibonding Cus-Pp close to the Fermi level. The larger footprint of Cus-Pp bonding contributions noticed in the –COHP for BLK-P in comparison to BL-P, compounded with the larger number of P-Cu bonds described earlier, suggests that these are the bonds that favor the adsorption of BLK-P over BL-P on Cu(311). It is noticeable that in case of BLK-P the hybridization of P${}_{p_{x}}$ orbitals is minimal and the bonding is dominated by P${}_{p_{y,z}}$ hybridization with Cu orbitals. This is expected from the orientation of P${}_{p_{y}}$ orbitals (perpendicular to the step), in agreement with Ref. Hoffmann (1988). In comparison, the BL-P $C_{3v}$ symmetry gives Pp orbitals that are neither strictly orthogonal or parallel to the step, resulting in a non-vanishing Cus-P${}_{p_{x}}$ contribution.

In conclusion, calculations within the DFT framework with vdW corrections suggest that BLK-P and BL-P are selectively formed on Cu substrates, depending on the choice of low- and high-index Cu surfaces. The high values of BLK-P and BL-P adsorption energies indicate strong interaction with the Cu substrates. The factors that favor BLK-P on stepped surfaces are the $C_{2v}$ symmetry of the surface and the increased DOS at the metal surface step. Here we used a simplified model system with a minimal terrace length, but we expect this to be valid for other metals and on wider terraces, albeit with different strain values. This work does not include alignment at different in-plane angles that may lead to Moire patterns. While the latter can lead to small variations in the adsorption energies on hexagonal and low-index surfaces, we expect that the $C_{2}v$ symmetry of BLK-P give significantly smaller adsorption energies for BLK-P unaligned with steps on large-index surfaces.

This work was supported in part by the DOD High Performance Computing Modernization Program and funded in part by the Air Force Office of Scientific Research under grant number FA9550-19RYCOR050. DH has been supported by an AFOSR-NRC postdoctoral fellowship through grant number FA9550-17-D-0001.