Strain-induced Weyl and Dirac states and direct-indirect gap transitions in group-V materials

We perform comprehensive density-functional theory calculations on strained two-dimensional phosphorus (P), arsenic (As) and antimony (Sb) in the monolayer, bilayer, and bulk $\alpha$-phase, from which we compute the key mechanical and electronic properties of these materials. Specifically, we compute their electronic band structures, band gaps, and charge-carrier effective masses, and identify the qualitative electronic and structural transitions that may occur. Moreover, we compute the elastic properties such as the Young's modulus $Y$; shear modulus $G$; bulk modulus $\mathcal{B}$; and Poisson ratio $\nu$ and present their isotropic averages of as well as their dependence on the in-plane orientation, for which the relevant expressions are derived. We predict strain-induced Dirac states in the monolayers of As and Sb and the bilayers of P, As, and Sb, as well as the possible existence of Weyl states in the bulk phases of P and As. These phases are predicted to support charge velocities up to $10^6$~$\textrm{ms}^{-1}$ and, in some highly anisotropic cases, permit one-dimensional ballistic conductivity in the puckered direction. We also predict numerous band gap transitions for moderate in-plane stresses. Our results contribute to the mounting evidence for the utility of these materials, made possible by their broad range in tuneable properties, and facilitate the directed exploration of their potential application in next-generation electronics.

In this work, our aim is to provide a comprehensive analysis of the monolayer, bilayer and bulk phases of orthorhombic P, As and Sb, in order to identify and compare the qualitative strain-related properties of each structure from a consistent set of calculations, thus treating each material on the same footing. Our findings provide new insights into their electromechanical properties, especially regarding arsenic and antimony, which have been relatively unexplored to date. Specifically, we identify qualitative transitions in band gaps, effective masses, structure, and topology that occur at various strains, and compute the elastic properties that determine the required stresses to attain these electronic states.
We predict the existence of strain-induced Dirac states in monolayer As and Sb, bilayer P, As and Sb, as well as possible Weyl states in bulk P and As, at moderate stress values. Our findings show that all of the predicted Dirac and Weyl points are indeed linear, at least in the Γ − Y direction, i.e. the puckered direction. Thus, following the convention of terminology found in the [38,48,71], and other sources, we classify Dirac or Weyl states as those associated with regions of sustained linear dispersion in the band structure, at or near the Fermi level, in at least one direction. We predict these states to support ballistic conduction and are unaffected by the spin-orbit coupling (SOC). In particular, few-layer P and As exhibit a strong indication of anisotropic conduction, dominated by ballistic conductivity along Γ − Y.
The outline of this article is as follows: in section 2 we review the details of our calculations and discuss the Voigt-Reuss-Hill averaging scheme used to compute the experimentally-relevant elastic properties. In section 3 we present the results of our calculations including the lattice constants, strain dependence of electronic properties, in particular the potential Dirac and Weyl states identified, and computed elastic properties. We conclude in section 4 with a discussion of the implications of our key results.

Calculation details
The calculations were performed with the QuantumEspresso package [72] using the Perdew-Burke-Ernzerhof (PBE) form of the generalizedgradient approximation (GGA) exchange-correlation functional [73]. An ultrasoft pseudopotential [74] from the SSSP Library [75] (with 5 valence electrons) was used to represent the core electrons. Non-SOC calculations were initially performed and those that exhibited potential Dirac or Weyl states were reassessed including non-perturbative SOC. In all calculations, van der Waals (vdW) interactions were incorporated using the B97-D empirical dispersion correction functional [76]. In order to achieve an energy convergence of at least 1 meV/atom and force convergence of at least × − 1.3 10 4 eV/a 0 , we found it sufficient to use a common plane-wave energy cutoff of 1100 eV with 'cold' smearing [77] of − 10 4 K for all elements. To achieve the same convergence, the Brillouin zone sampling for bulk systems was × × 15 15 15, and × × 15 15 1 for monolayers and bilayers. Uniaxial and shear strains between ±5% were applied in increments of 1% to the unit cell with internal relaxation subject to the same force convergence criterion as above. Electronic band structures were calculated along the high-symmetry points of the Brillouin zone {Γ, X, S, Y, Z} (figure 1(b)) for each value of in-plane strain (figures S1−S9 in the supplemental material (stacks.iop.org/TDM/4/045018/mmedia)). For shear strains the Brillouin zone deforms into an asymmetric honeycomb, yet we continued to sample along the original path since the deformation up to 5% strain is negligible and the effective masses are all  calculated at the Γ-point. We determine the Kohn-Sham band gap from the band structures and chargecarrier effective masses according to the nearly-free electron model 1 using a cubic spline fit to 9 data points about the Γ-point. The charge velocities were similarly determined according to the dispersion relation from the linear fit to the Dirac or Weyl lines. The elements of the stiffness matrix C were derived from the gradients of the resultant stress-strain profiles / σ ε = ∂ ∂ c ij i j , from which all elastic properties were derived. In practice, however, the computed stiffness tensors are not exactly symmetric due to numerical noise but we make them so by taking the average of C and its transpose C T as the effective stiffness tensor. We begin our discussion with a brief overview of the the Voigt-Reuss-Hill scheme, which is a popular model used for computing effective isotropic elastic properties.

The Voigt-Reuss-Hill scheme
In order to effectively preserve, study and strainengineer few-layer nano-structures, such as BP [78], graphene [79], or molybdenum disulfide [80] (MoS 2 ), the nano-flakes are typically deposited onto a suitable substrate. The cumulative contribution of dispersed nano-flakes distributed on or within a bulk medium results in the macroscopic elastic properties that are measured by experiments. The theoretical calculation of these elastic properties requires an appropriate mixture model (such as the rule-of-mixtures [11] (ROM) or the Halpin-Tsai [81] models) that require the (typically averaged) elastic properties of the interstitial nano-flakes.
In the theory of effective media, isotropic bulk properties are computed by averaging the stiffness tensor C over all possible rotated reference frames [82][83][84], as outlined in the supplemental material. The result is called the Voigt average [85,86] and it gives isotropic averages for the bulk Young's modulus Y V , and shear modulus G V , given in equation (S7). The same scheme applied to the compliance tensor S results in the corre sponding Reuss averages [87], Y R and G R , given in equation (S9). The Voigt scheme assumes that the material undergoes constant strain and it returns over-estimated elastic constants. Conversely, the Reuss scheme assumes constant stress and it tends to underestimate the elastic constants. The Hill averages [88] , and 2 , from which the isotropic Poisson's ratio ν H and bulk modulus B H are expressed as These are widely considered as reliable estimates of the actual physical values [84]. The Voigt-Reuss-Hill approach described above is used in the present work to determine the isotropic averages of the Young's, shear and bulk moduli, and the Poisson ratio of the bulk P, As, and Sb structures using the elements of the elastic tensors. Let us now discuss how the above approach may be adapted to derive the relevant equations for the specific case of two-dimensional materials.

In-plane Voigt-Reuss-Hill average
If the interstitial nano-flakes in a bulk medium form high-quality planar sediments [11,[89][90][91][92][93][94][95], the random orientation occurs instead in the plane of the flakes and we must calculate isotropic-averages in-plane. Due to the weak vdW bonds between layers, strains related to out-of-plane directions can be ignored resulting in the reduced-stiffness tensor (equation (S3)). In the supplemental material, we re-derive the angular dependence of the rotated tensor-elements ( ) θ C ij and ( ) θ S ij about the → z-axis (equation (S6)) as a function of the elements in the original reference frame, similar to the general Voigt-Reuss scheme. The angulardependence of the in-plane elastic constants are then expressed as (3) with the Hill-average taken as in equation (1). By integrating the elastic tensors ( ) θ C ij and ( ) θ S ij over π 2 , the in-plane averages are then computed analogously.

Lattice constants
The lattice constants a, b, c of the fully-relaxed structures are presented in table 1, where, in the Table 1. Lattice parameters (Å) for monolayer, bilayer and bulk structures of P, As, Sb compared to experimental data [100][101][102] quoted in parentheses. For the monolayers and bilayers the layers thickness ′ c is given. monolayer and bilayer cases, we quote the layer thickness ′ c instead of the unit cell height c. Our computed lattice parameters compare well with other recent theoretically predicted values [28,29,51,59,[96][97][98][99] and the available experimental data [100,101]. For a given element, we find that the lattice parameter along the puckered direction, 'a', shortens as the number of layers increases, which agrees with observations in other studies. This is attributed to the increased vdW forces between layers, which leads to increased binding primarily in the softer puckered direction.

Electronic properties
All of the band structures pertaining to the following analysis are presented in figures S1−S9 of the supplemental material. Where we identify possible Dirac or Weyl states, high resolution, three-dimensional band structures with SOC at representative strains are recalculated. To confirm the existence of linear-dispersion, we also plot lines along the surface of the Dirac and Weyl points at 0 , 30 , 60 and 90 with respect to the Γ − X line. A representative sample of these results are presented in figures 5(a)-(g), while the rest can be found in figures S10−S13 of the supplemental material.
In general we find the band gap to be very responsive to uniaxial in-plane strain but significantly less so with respect to shear strain. We identify several directindirect band gap transitions as well as the opening and closing of band gaps, summarized in table 3. We find the charge-carrier effective masses vary approximately linearly with respect to the uniaxial strain in general but with notable exceptions that will be discussed. This section is divided into three parts discussing each of the species-P, As and Sb-for which we review the qualitative calculation results including band gap trans itions, effective mass behavior and linearly-dispersive bands.

Phosphorus
As shown in figure 2(a), our calculations reproduce the direct band gap of 0.88 eV at the Γ-point in the relaxed P monolayer, which falls within range of the reported gap between 0.7 eV (DFT-PBEsol [12]) and 1.0 eV (DFT-HSE06 [8]). On the other hand, quasiparticle calculations predict a larger 2 eV band gap [22] with significant exciton binding [103] (between 0.4-0.83 eV). However, it is well understood that approximate semi-local functionals such as PBE suffer from a systematic band gap problem [104] that may also adversely affect the metal-insulator critical strains. Nevertheless, it is important to emphasize that band alignments and rates of change are quite often reliably reproduced [105], as are the direct-indirect transitions [13] in two-dimensional materials. While absolute band gaps are therefore not expected to be exactly reproduced, we can expect reasonable agreement with trends in electronic and mechanical behavior [106,107]. The application of uniaxial in-plane strain is found to open the band gap for tensile strain and diminish it for compressive one, while shear strain has a negligible effect. The electron and hole effective masses (figures 2(b) and (c)), compare well to the figures computed in [108], where, at ε = +5% xx tensile strain, the electron and hole effective masses coincide at 0.9 m 0 as higher energy bands fall below the conduction band. For compressive strains, the hole effective mass along Γ − X rises significantly as the valence band flattens.
Bilayer P is also found to have a direct band gap of 0.43 eV (figure 2(d)) in the relaxed state, and broadly the same behavior as the monolayer, in which case the band gap closes at around −5% uniaxial compressive strain with a predicted Dirac state at the Γ-point.  [54], who demonstrated the effect of strain on hopping parameters can lead to a Dirac semi-metallic state in bilayer P. Similarly, Baik et al [71] found that the SOC did not induce a band gap in potassium-doped multi-layer P, but did lift the spin-degeneracy of the Dirac points. A directindirect band gap transition is also observed at +2% uniaxial tensile strain. The effective masses (figures 2(e) and (f)), also exhibit broadly the same behavior as the monolayer, due to band-flattening at Γ − X and the falling conduction bands along Γ − Y, which lead to the charge carrier effective masses along Γ − X converging at +4% strain and an increasing hole effective mass for compressive strains.
In the bulk, however, we find that the band gap is completely closed (figure 2(g)), i.e. that the material is metallic. After investigation, we concluded that this was an effect of the smearing functionality [77] in the relaxation procedure and that it contradicts numerous experiments [109][110][111][112] that have measured a direct gap in the range of 0.31-0.36 eV. When relaxed under fixed-occupancy conditions, instead, a band gap of ∼0.35 eV was produced. While the PBE gap remains closed in the relaxed state, under uniaxial tensile strain it briefly becomes a single-point semi-metal at +2%. At such strains a possible Weyl state is observed, before a direct gap opens that subsequently transitions to an indirect gap at +3%. Shown in figure S10 in the supplemental material is the three-dimensional band structure with SOC in which a pair of potential Weyl points occur on an off-symmetry point ′ X along Γ − X. Here, the SOC slightly reduces the band gap by ∼0.05 eV and does not qualitatively affect the overall results. The maximum charge velocity is ( ) = × v 2.40 1 10 6 m − s 1 for both ε = −5% xx and ε = −5% yy and occurs along a line parallel to Γ − Y. Under greater compression this band-inversion may also lead to further Weyl states, which have been experimentally observed at similar pressures [60][61][62]. Again, shear strain is seen to have a negligible effect on the gap. The electron effective masses are quite responsive to strain (figure 2(h)), where those along Γ − Y rise for both tensile strain along ε xx , due to falling conduction bands, and compressive strain along ε yy due to flattening bands along Γ − Y. The effective masses along Γ − X were necessarily not computed once the bands overlapped below +2% strain and the hole effective masses are found to vary with respect to the strain to a slightly lesser extent ( figure 2(i)).
To summarize, we predict the onset of Γ-point Dirac states in bilayer P at −5% uniaxial compressive strain, with effective one-dimensional conductivity at ε = −5% xx , and a direct-indirect band gap transition at +2% tensile strain. We also predict the existence of a possible Weyl states at +2% tensile strain in bulk P, followed by a direct-indirect band gap transition at +3%. Finally, effective masses are found to be particularly responsive to ε xx uniaxial strain.

Arsenic
In contrast to P, we identify an indirect band gap of 0.15 eV along the Γ − Y direction in the relaxed As monolayer ( figure 3(a)), which is significantly lower than the predicted DFT-HSE06 gap [98] of 0.83 eV. However, the relaxed band structure and band gap profiles closely resemble those in [28,49]. The band gap diminishes for tensile strain along ε xx and at +2% the material becomes semi-metallic with a Dirac state at the Γ-point emerging at ε = +5% xx accompanied by an electron pocket above the Fermilevel (figure S11 in the supplemental material), which is unaffected by the SOC. The maximum charge velocity here is ( ) = × v 3.01 1 10 6 m − s 1 and lies along Γ − Y . In the orthogonal direction the bands are flat, similarly to monolayer P, with a relatively small charge velocity. This high anisotropy in charge velocities, dominated by the ballistic conduction along Γ − Y , is again indicative of effective one- For compressive strain along ε xx an indirectdirect transition occurs [28] at ε = −2% xx . For tensile strain along ε yy the band gap opens, where at ε = −3% yy , the indirect band gap closes along Γ − Y. Similar to monolayer phosphorus, there is no appreciable effect due to shear-strain. Meanwhile, the chargecarrier effective masses (figures 3(b) and (c)) respond linearly to uniaxial strain and compare well to other works [101], where, in particular, valence band broadening along Γ − X leads to an increasing hole effective mass.
For bilayer As, we identify a direct band gap of 0.45 eV ( figure 3(d)), in contrast to the indirect band gap observed in the monolayer. Here, the band gap opens for uniaxial tensile strain and diminishes for compressive strain. The direct band gap transitions to an indirect gap at both ε = −3% yy and ε = +2% yy , while at ε = +2% xx it also transitions to an indirect gap before resuming to a direct gap again at ε = +3% xx . Moreover, we predict a Dirac state at the Γ-point at a compressive strain of ε = −4% xx (figure S11 in the supplemental material) for which the maximum charge velocity is ( ) = × v 2.62 2 10 6 m − s 1 along Γ − Y. Along Γ − X the bands are also flat, similar to the monolayer, and have a relatively negligible charge velocity. The high anisotropy in charge velocities, is again indicative of effective one-dimensional conductivity, dominated by the ballistic conduction along Γ − Y, and is further supported by the large disparity in effective masses at ε = 5% xx , shown in figures 3(e) and (f). Here again, the SOC has no appreciable effect on the bands. The electron and hole effective masses (figures 3(e) and (f)) respond approximately linearly to the applied strain, where conduction band broadening leads to increased effective electron masses, and valence band flattening at Γ leads to increasing hole effective masses for strain along ε yy .
Finally, no band gap is determined in the relaxed bulk phase (figure 3(g)), again contrary to experiments [113], where a small direct band gap of ∼0.3 eV is observed. However, at ε = +1% xx strain, a potential Weyl state is briefly observed on an off-symmetry point ′ X (figure S12 in the supplemental material) before a direct gap opens that subsequently trans itions to an indirect one at ε = +3% xx , after which it reduces again. Another potential Weyl state around the same off-symmetry point ′ X is also predicted to occur between ⩽ ⩽ ε + + 1% 2% yy after which a direct band gap also appears. The recalculated band structure with the SOC for ε = +1% yy (figure S12) confirms the linear-dispersion. The maximum charge velocity in  In summary, we predict Γ-point Dirac states in the monolayer and bilayer of As, which support onedimensional ballistic conduction, as well as possible Weyl states on off-symmetry points in the bulk at moderate levels of in-plane stress that are unaffected by the SOC. We also observe several band gap transitions, in particular in the monolayer phase, which also include semi-conducting-metallic transitions. Finally, the effective masses respond approximately linearly with respect to uniaxial strain, except in the bulk, which exhibits quadratic behavior.

Antimony
The relaxed Sb monolayer is found to possess an indirect band gap of 0.21 eV along Γ − Y ( figure 4(a)), which is reasonably comparable to other PBE values 0.28 [29]-0.37 [114] eV, although these have been obtained by including SOC. For tensile strain along ε yy the band gap opens, suggesting an indirect-direct transition for strains above − 6 7%, and it diminishes for compressive strains before finally closing at ε = −2% yy , where the material becomes a semi-metal. Similarly, the indirect gap closes along Γ − X at a compressive strain of ε = −2% xx at which monolayer Sb again becomes semi-metallic. The indirect gap transitions to a direct gap at ε = +1% xx tensile strain and remains so until finally closing at ε = +4% xx , at which point we predict a potential Dirac state along Γ − Y at an off-symmetry point [37] ′ Y ( figure 5(g)) that has also been predicted in [51].  5(i)). Moreover, the valence band at the X-point undergoes a Rashba splitting [115] due to SOC, which is also predicted to occur in the monolayers of α-P [116], and β-Sb [52,117]. Finally, the electron effective masses experience a weak linear response to strain (figures 4(b) and (c)), while the hole effective masses along Γ − X respond much more strongly to a rapid broadening or flattening of the valence band.
Furthermore, the relaxed bilayer phase is found to be semi-metallic where an indirect band gap opens at ε = +3% yy tensile strain and for uniaxial strains < −1% band-inversion at the Γ-point leads to to a fully-metallic state. In addition, a possible Dirac state emerges at a similar non-symmetry-point ′ Y along Γ − Y for ε = +2% xx tensile strain (figure S13 in the supplemental material) and remains in place up to at least +5% strain. The maximum charge velocity ( ) = × v 4.47 3 10 6 m − s 1 is also along Γ − Y and is approximately the same as that of the monolayer. The effective masses (figures 4(e) and (f)) experience mild linear-response to strains prior to the transition to full metallicity, at which point a rapid flattening of the bands at the Γ-point suggesting strong electron localization Beyond a compressive strain of ε = −3% yy , however, at a stress of ∼0.3 GPa, the bilayer undergoes a structural transition and buckles in the puckered ( → y) direction. This buckled structure has a total energy 1.7 meV/atom lower than that of the relaxed state of the unperturbed α-bilayer and 3.0 meV/atom lower when allowed to fully-relax, as shown above in figure 1(c). This suggests the possible existence of a new structure that is attainable via strain. Finally, bulk Sb is found to be completely metallic for all levels of strain explored in this work. However, shear strains in this case do appear to have a significant effect on the bands despite not opening a gap. In summary, we predict possible non-symmetry-point Dirac states in the strained monolayer and bilayer of Sb, which are qualitatively unaffected by SOC, as well as Rashba splitting at the X-point in the monolayer. We also predict indirect-direct and indirect-semi-metallic transitions in the monolayer phase and a band gap opening in the bilayer phase. Finally, we observe a buckled state induced in bilayer Sb at −4% compressive strain. Bulk Sb was found to be metallic at all levels of strain explored.
In table 2 below we present a summary of the calculated band gaps and effective charge carrier masses for the relaxed phases of each structure, and in table 3 we provide a synopsis of the band gap and phase transitions of interest. Finally, in figure 5, we present the band structures of bilayer P and monolayer Sb, that form a representative sample of the different Dirac points predicted at Γ, ′ X and ′ Y , as well the three-dimensional band structures about the region of the points.

Isotropic bulk properties
In order to obtain these electronic states, and to ensure accurate strain-engineering, knowledge of the mechanical properties is paramount. Therefore, in this section we review the mechanical response of the few-layer and bulk phases in order to compute the elastic properties, both isotropically averaged and as a function of orientation of applied in-plane stress.
The computed elements of the stiffness tensor C in GPa of each structure are given in table SI in the supplemental material, where those pertaining to bulk P compare well to experiments [118,119] and similarly computed values [96,108]. For the elements related to in-plane strains (c 11 , c 22 , c 66 , c 12 ) we observe an expected increase in stiffness as the layer number increases and for decreasing atomic number. However, for other elements relating to out-of-plane and shear stresses (c 33 , c 44 , c 55 , c 23 , c 13 ) the stiffness actually increases in the bulk phase from P to As to Sb.
The Hill-averaged bulk properties are presented in table 4, which compare well to other DFT values [96,108], while our calculated bulk modulus for bulk P (37.2 GPa) is also within reasonable range of the experimental values (32.32 [100]-36.02 GPa [120]). We also observe that the bulk properties remain largely comparable for all the species, but generally decrease from P to As to Sb (except for the Poisson's ratio and bulk modulus, which are largest for As). We also note that while bulk P has the largest in-plane responses, As and Sb have larger out-of-plane and shear responses, which enable the net isotropic properties for all the species to remain comparable overall.

In-plane elastic properties
In section SI of the supplemental material we re-derive the equations for the elastic properties as a function of the in-plane orientation angle θ, defined in figure 1, as outlined in [121]. These functions are the plotted in figure 6 and include the Young's modulus ( ) θ Y and it's average ⟨ ( )⟩ θ Y , the shear modulus ( ) θ G , and Poisson's ratio ( ) ν θ . The experimental Young's modulus of 130 GPa, determined in [11] via the ROM of a nanoflake-polymer composite, lies precisely between the average in-plane bulk Young's modulus 85.7 GPa and the elastic stiffness in the zigzag ( → x) direction 187.9 GPa, thus fitting the results of our model reasonably well.
The anisotropy of the elastic properties is also apparent in the mechanical profiles, particularly with regard to the Young's modulus, which has a 2-fold symmetry about the x-axis, in contrast to the shear modulus and Poisson's ratios, which display 4-fold symmetry about both the axes (except for the Poisson's ratio of P, which remains 2-fold symmetric). Indeed, for a given species, the general shape of each profile is approximately preserved with respect to the number of layers, while the range of each property tends to increase. This discovery is advantageous in the strain-engineering of nano-flakes since one may now forecast in advance the response of a material to inplane strain, once the underlying profile and number of layers are known.
Another interesting feature is that the extrema of the elastic functions do not necessarily coincide with the coordinate-axes. For instance the Young's modulus maximum for monolayer Sb occurs at 22 . Table 5 summarizes the global minima and maxima of each function and the angles at which they occur. While most of the function extrema occur expectedly at 0 , 45 or 90 , many are incident away from the coordinate-axes. This result lends further insight into the Table 2. Kohn-Sham band gaps (eV), indicating the ( a ), ( b ) and ( c ) states, as well as the charge-carrier effective masses (m 0 ) for each phase of P, As, Sb.
mechanical aniso tropy of the orthorhombic group-V materials.
The emergence of a negative Young's modulus of −2.1 GPa in monolayer As ( figure 6(b)), at first glance, may give cause for concern. It arises due to a negative Voigt estimate for the Young's modulus at 90  GPa (since > c c 12 22 ), which is larger in absolute magnitude than the Reuss estimate at 90 given by / = s 1 4.7 22 GPa and results in a net negative Hill-average. In this instance, we surmise Table 4. The Hill-averaged Young's modulus Y H , shear modulus G H and bulk modulus B in GPa, and Poisson's ratio ν H for bulk P, As and Sb compared to similarly calculated DFT [96] values and available experimental data [100,120].  that either the assumptions of the Voigt model break down, or the Hill-method is not universally appropriate in arbitrary directions in-plane and a more robust averaging scheme must be employed. Nevertheless, the qualitative in-plane functions and their isotropic averages remain physically meaningful and in general can provide valuable physical insight.
In contrast to the isotropic averages for the bulk properties in table 4, which remained largely comparable for P, As and Sb, a much clearer trend across the species emerges once we have eliminated contributions from the out-of-plane and shear stresses. In this instance, P clearly possesses superior in-plane mechanical strength in both moduli, which decrease with the number of layers as expected As and Sb are largely similar in the monolayer, though less so in the bilayer and bulk phases, where they are stronger in the → x-direction. In contrast, the Poisson's ratio tends to remain relatively stable aside from generally decreasing with increasing number of layers.
In summary, there exists in the elastic properties a broad range of responses, profile shapes and behavior that is reflective of the underlying anisotropic crystal structure, which also strongly depend on the number of layers with P typically the stiffest and Sb the most flexible. In general we find the shapes of the in-plane response profiles to be conveniently consistent, which implies that, for a given number of layers, the in-plane elastic response of a nano-flake can be reliably estimated a priori.

Conclusion
We have extensively explored the mechanical and electronic properties of P, As and Sb in their few-layer and bulk phases. We have identified several band gap transitions in almost all of the structures. The SOC tended to close the bands by ∼0.05 eV but did not alter any of our qualitative findings. We also predict the existence of Dirac states in the strained phases of monolayer As and Sb, bilayer P, As and Sb as well as possible Weyl states in bulk P and As for moderate levels of strain. The linear-dispersion was observed along Γ − Y of each of the predicted Dirac or Weyl states, corresponding to the direction of softest mechanical response in the puckered direction. The maximum charge velocity is calculated to be over 10 6 m − s 1 . In particular, for bilayer P and few-layer As we predict highly anisotropic conductivity dominated by ballistic transport along the puckered direction that is indicative of effective, one-dimensional conduction. We predict that an appropriate strain could yield these effects in experiments.
We also observe the existence of a notable buckled state of compressed bilayer Sb at −4% strain. Finally, the angular-resolved elastic properties as well as the stress-dependence of the Kohn-Shame band gaps and charge-carrier effective masses revealed highly anisotropic behavior, spanning a broad range of values, and angular-dependent behavior that has become characteristic of these group-V layered structures. Moreover, the critical stresses at which these transitions occur are expected to be experimentally accessible and highly switchable, paving the way for possible verification in the near future. Thus, the group-V layered materials are poised to become central to the next generation of electronic devices with potential novel applications in field-effect transistors; batteries; gas-sensors and opto-electronic devices. the work conducted for this Article. This work was enabled by Science Foundation Ireland (SFI) funded centre AMBER (SFI/12/RC/2278). All calculations were performed on the Kelvin cluster maintained by Trinity College Dublin Research IT and funded through grants from SFI.