The freestanding monolayer BP has a graphene-like planar structure due to the sp2 hybridization of atoms. After optimization, the B–P bond in monolayer BP has a length of 1.856 Å, which is longer than the B–N bond length of 1.45 Å in monolayer BN [37], because P atom has a larger atomic radius with respect to the N atom. The calculated in-plane lattice constant is 3.215 Å, which is in agreement with previous theoretical calculations [12, 37, 38].

For bilayer BP, the most stable configuration is shown in Fig. 1. The B atoms on top layer are located above the P atoms on bottom layer, while P atoms on top layer centered above the hexagon of the bottom layer. The B-P bond length for bilayer structure does not change with respect to the monolayer structure. The optimized interlayer distance is 3.451 Å which is in agreement with previous studies [39, 40].

Figures 2 (a) and (b) display the calculated electronic band structures of monolayer and bilayer BP along the high symmetry points of the first Brillouin zone. The monolayer BP has a direct band gap and both the Conduction Band Minimum (CBM) and Valence Band Maximum (VBM) locate at the K point. For bilayer BP, the CBM is located at the K point while the VBM is slightly shifted from the K point and it has an indirect band gap. The calculated band gap of monolayer and bilayer BP is 0.90 eV and 0.72 eV, respectively. Not that the band gap of bilayer structure is smaller than the bandgap of monolayer structure because of the inter layer interactions. The smaller band gap of bilayer BP compared with that of the monolayer BP makes it more suitable for thermoelectric applications. The calculated partial density of states for monolayer and bilayer BP are illustrated in Figs. 2 (c) and (d). It can be seen that for both structures the VBM is formed mostly by the p orbital of P atoms. For monolayer BP the CBM is mostly composed of p orbital of B atoms while for bilayer BP the CBM is equally distributed with the p orbitals of B and P atoms.

Figure 3 (a) shows the phonon dispersion of monolayer BP. As can be seen from Fig. 3 (a), there is no imaginary frequency, which shows the dynamic stability of monolayer BP. It can be seen that as **k** *→* 0 the phonon dispersions of LA and TA modes are linear but that of ZA mode is quadratic. The separation of optical LO and TO modes from three acoustical modes shows an indirect frequency gap of about 33 meV.

To study the effect of interlayer interaction, we also plot the phonon dispersion of the bilayer BP in Fig. 3 (b). There are 12 phonon modes in the bilayer structure, while there are 6 modes in the monolayer structure.

If two layers are placed together without interaction, we expect that their phonon dispersions completely overlap. The phonon dispersion of bilayer is very similar to the one of monolayer BP and most of the branches in the bilayer can be viewed as the ones split from the monolayer. Almost no change is observed in the in-plane doubly degenerated bands (the in-plane TA, TO, LA and LO) while, the interlayer coupling causes the splitting of out-of-plane ZA and ZO modes at the Γ point. The ZA and ZO1 modes show more obvious splitting with respect to the ZO2 and ZO3 modes.

Next, we study the phonon dispersion properties of BP sheet under in-plane biaxial strain. As we can see in Fig. 4, the ZA phonon mode becomes imaginary and softens near the Γ point when the applied compressive strain is − 2%. By increasing the compressive strain, the larger absolute values of the maximum imaginary frequency becomes obvious. Therefore, applying the compressive strain over − 2% leads to dynamic instability. Also, the LO/TO modes at Γ-point becomes harden under compressive strain. The phonon bandgap between TO/LO phonon modes with the LA phonon mode at Γ-point increases, as compressive strain increases. As shown in Figs. 4 (f)-(j), the monolayer BP is stable under tensile strain because no imaginary frequency is found under the tensile strain. It can be observed that the ZA phonon mode near Γ-point changes from quadratic to linear. We also notice that under biaxial tensile strain, all phonon modes except of ZA mode away from Γ-point become soften. This is because of the weakness of atomic bond strength due to the increasing of lattice constant by applying the tensile strain. Moreover, by increasing of tensile strain, the LO/TO modes at Γ-point becomes soften, which results in the reduction of phonon bandgap between LA and LO/TO modes.

Figure 5 shows the phonon dispersion of bilayer BP under biaxial tensile strain. When a tensile strain is applied to bilayer BP, the phonon modes, specifically LO2, LO3/TO2, and TO3, experience a shift towards lower mode frequencies. This shift results in phonon softening and a reduction in the phonon band gap between (LO2, LO3/TO2, TO3) and LA and LO1. Additionally, the ZA mode undergoes a hardening process, which is similar to what occurs in the monolayer structure. Overall, these changes in the phonon modes can significantly impact the thermal properties of bilayer BP.

**Optical properties**

After examining the phonon dispersion properties of bilayer BP and monolayer BP without and with strain, the frequency dependent real ε1(ω) and imaginary part ε2(ω) of dielectric function, absorption coefficient α(ω), and energy loss function L(ω) for all structures are investigated and demonstrated in Figs. 6–11.

For x direction of polarization, the unstrained monolayer BP has two absorption peaks at 2.65 eV and 6.64 eV and the bilayer BP has two peaks around 2.95 eV and 6.76 eV, respectively (see Figs. 6(a) and 7 (a)). The onset of absorption is located around ~ 0.8 eV and ~ 0.6 eV for monolayer and bilayer BP, respectively, which is in compliance with a direct band gap at the K point with energies, 0.90 eV and 0.72 eV, respectively. For z direction of polarization, the unstrained monolayer BP (bilayer BP) has main absorption peaks at 8.42 eV, 9.95 eV, 12.61 eV and 14.84 eV (2.90 eV, 7.74 eV, 9.78 eV and 14.11 eV). The calculated absorption coefficient of monolayer BP indicates that the monolayer BP has high sensitivity in the ultraviolet regime. For x direction of polarization, the ε2(ω) for monolayer BP (bilayer BP) shows significant peaks at 1.16 eV, 2.52 eV and 6.52 eV (1.37 eV, 2.19 eV and 6.52 eV). The real part, ε1(ω), denotes the electronic polarizability of the material and dispersion effects. Figures 8 (c) and (d) and 9 (c) and (d) show ε1(ω) for light polarized along the x and z directions for monolayer and bilayer BP, respectively. The magnitude of ε1(ω) at ω = 0 is the most important quantity of this part of dielectric function. Without the applied strain, the ε1(0) value of monolayer BP is calculated to be 4.03 and 1.49, along the x and z directions. The calculated ε1(0) for bilayer BP, is found to be 7.66 and 2.47, along the x and z directions.

For x direction of polarization, the absorption edge shifts slightly to lower energies and the main absorption peaks of monolayer BP become blue-shifted and their intensities increases and the value of ε1(0) increases as the strain changes from ε = +10% to ε = -10%. Also, the first peak of ε2(ω) becomes red-shifted, while the second and third peaks become blue shifted as the strain changes from ε = +10% to ε = -10%. Therefore, application of compressive strain increases the absorption of BP monolayer by increasing the intensity and width of absorption.

The electron loss function for monolayer and bilayer BP is demonstrated in Figs. 10 and 11. For x direction of polarization, the main peaks of L(ω) for monolayer BP were detected at 3.33 eV and 9.56 eV. The first one spreads from the infrared to the near ultraviolet region. The second band is found mainly in the middle and far ultraviolet region. In the electron loss function of bilayer BP, we have found dominant peaks at 3.61 eV and 12.23 eV. For monolayer and bilayer BP, it can be predicted that the main resonance and electron energy loss occurs in the ultraviolet regime. For z direction of polarization, the main peaks of L(ω) for monolayer BP (bilayer BP) were observed at 8.94 eV, 10.46 eV, 12.74 eV and 15.20 eV (3.0 eV, 8.34 eV, 11.01 eV, 14.79 eV and 16.01 eV). It can be observed that the number of peaks in the bilayer structure is more than the monolayer structure, due to the interaction between the layers along the z direction.

Figure 10 displays the electron loss function of monolayer BP for various biaxial strains. For x direction of polarization, it can be seen that by increasing the compressive (tensile) strain the position of main peaks is blue (red) shifted and the intensity of first main peak increases (decreases), while the intensity of second peak decreases as the changes from ε = +10% to ε = -10%. For z direction of polarization, the position of all peaks is blue shifted when the strain varies from ε = +10% to ε = -10%. Figures 7, 9 and 11 show that the general behavior of absorption coefficient, real and imaginary part of dielectric function and energy loss function for bilayer BP under small strain is similar to monolayer BP. For x direction of polarization, the main peaks of α(ω) and L(ω) are blue shifted and their intensities increases as the strain changes from ε = +10% to ε = -10%. Also, by applying the strain from ε = -10% to ε = +10%, the distance between the dual neighbor peaks in the low energy region of ε2(ω) decreases and the main peak in the high energy region is red shifted.