The diverse magneto-optical selection rules in bilayer black phosphorus.

The magneto-optical properties of bilayer phosphorene is investigated by the generalized tight-binding model and the gradient approximation. The vertical inter-Landau-level transitions, being sensitive to the polarization directions, are mainly determined by the spatial symmetries of sub-envelope functions on the distinct sublattices. The anisotropic excitations strongly depend on the electric and magnetic fields. A uniform perpendicular electric field could greatly diversify the selection rule, frequency, intensity, number and form of symmetric absorption peaks. Specifically, the unusual magneto-optical properties appear beyond the critical field as a result of two subgroups of Landau levels with the main and side modes. The rich and unique magnetoabsorption spectra arise from the very close relations among the geometric structures, multiple intralayer and interlayer hopping integrals and composite external fields.

Introduction

Phosphorene, a two-dimensional allotrope of the group-V element black phosphorus (BP), has recently attracted many theoretical and experimental studies. Layered phosphorene presents the puckered configuration, which is different to the planar structure in graphene. It is closely related to the sp3 hybridization of four (3s, 3p, 3p, 3p) orbitals in BP. Furthermore, the deformed hexagonal lattice in the x-y plane is quite different from the honeycomb lattice of group-IV systems[1]. This unique geometric structure fully dominates the low-lying energy bands which are highly anisotropic in the dispersion relations as a function of the in-plane wave vector k, i.e., the linear and parabolic dispersions near the Fermi energy E, respectively, along the and directions[1]. The carrier effective masses are anisotropic and asymmetric between electrons and holes[2,3]. Accordingly, the best functional performances in nanoelectronics and optoelectronics depend on the rotation-induced orientations of BP systems, which could be deduced from optical spectroscopy measurements[3]. The strong anisotropy sets phosphorene aside from most 2D materials, such as graphene, boron nitride and transition metal dichalcogenides.

Up to now, phosphorene has been successfully synthesized by various experimental techniques, covering mechanical cleavage[4-6], liquid exfoliation[7-10], mineralizer-assisted short-way transport reaction[11-13], growing on red phosphorus (RP)[14], converting from RP[15], pulsed laser deposition[16] and plasma thinning[17]. A layer-dependent energy gap of ~0.5–2 eV[1,18,19] in phosphorene has been verified by optical experiments[5,20]. Such band gaps are larger than that of its bulk counterpart (~0.3 eV)[4,18,21] and are in sharp contrast with the zero or narrow gaps of two-dimensional (2D) group-IV materials[22]. Transport measurements have demonstrated that the phosphorene-based field-effect transistor exhibits an on/off ratio of 105 and the carrier mobility at room temperature as high as 103 cm2/V⋅s[4]. The high performance implies potential applications in nanoelectronics. Specifically, many remarkable anisotropic physical phenomena are reported, such as mechanical strains[23], excitonic effects[2], optical spectra[4]; charge transport[24-26] and thermal[27] properties. The peculiar characteristics can be traced back to its feature-rich lattice structure and electronic properties.

Phosphorene exhibits unusual optical properties, such as high optical absorption in the UV region and dichroism. The threshold absorption structure, corresponding to the band-gap energy, is revealed with light polarization along the armchair direction[3]. Its frequency falls off rapidly with the number of layers. In contrast, the band-gap absorption is forbidden under the zigzag-direction electric polarization, mainly owing to the specific symmetry of wave functions in the initial and final states. Theoretically speaking, the electric dipole excitations, which connect the valence and conduction band states, are completely different for two perpendicular polarization directions[3]. The strong anisotropy is also predicted in the magneto-optical conductivity of BP thin films[28]. Interestingly, the essential physical properties of BP could be greatly diversified by a uniform electric field (E)[29,30]. E might create a gapless band structure after reaching a critical value (E) or induce various energy bands such as parabolic bands, graphene-like Dirac cones and oscillatory bands. This arises from a competitive or cooperative relation between the intralayer and interlayer atomic interactions and the Coulomb potential energy. The strong electrically tunable energy bands enrich the quantization phenomena under a uniform perpendicular magnetic field (B), including two subgroups of Landau levels (LLs), uniform and non-uniform LL energy spacings and frequent crossings and anti-crossings[31]. There exist dramatic changes in the sub-envelope functions for the two mixed LLs with anti-crossing behaviors.

Motivated by the above-mentioned works, here we investigate the tuning effects, brought about by a perpendicular electric field, on the magneto-optical absorption spectra of bilayer phosphorene. The relations among the geometric structure, intrinsic interactions and external fields are explored in detail. The generalized tight-binding model[32] has been developed to investigate the electronic properties of monolayer and bilayer phosphorus in a composite electric and magnetic field. The electric-dipole transition elements are evaluated within the gradient approximation[33,34]. There exists certain challenges in performing the numerical calculations of the magneto-absorption spectra of bilayer phosphorene. Since the magnetic Hamiltonian is a very huge matrix within achievable experimental field strengths (details in the Methods section), a long period of time is needed for solving the eigenvalues and eigenvectors. Moreover, the matrix elements of the velocity operators (Eq. (3) in the Methods section) are hard to compute. To overcome these challenges, we utilize a band-like matrix[32] and the strong localization features of Landau wavefunctions to efficiently derive the absorption spectra. In addition, we use the method of gradient approximation to evaluate the matrix elements of the velocity operators as shown in Eq. (4).

In order to thoroughly explore the magneto-optical properties which significantly reflect the main features of LLs, we first present a simple review of magneto-electronic properties for which we refer to Figs 1 and 2[31]. Moreover, the effects of an electric field are demonstrated by the E-dependent LL wavefunctions (Fig. 2(b–g)). The magneto-optical absorption spectra are predicted to be strongly dependent on the polarization directions and the electric- and magnetic-field strengths (demonstrated in Figs 3–7). Through careful calculations and detailed analysis, the available inter-LL transitions are deduced to be mainly determined by the spatial distributions/symmetries of sub-envelope functions on the distinct sublattices. The selection rule, frequency, intensity, number and structure of absorption peaks are dramatically changed by the E field. Whether the usual/unusual optical properties are revealed in the magneto-absorption spectra relies on the electric field beyond the critical one or not. This is well explained by exploring the diverse selection rules. These unique phenomena, occurred in the near to mid-infrared region, might be important for a wide range of optical technologies, such as spectroscopy, materials processing and chemical and biomedical sensing. They result from the peculiar geometric structures, multiple intralayer and interlayer hopping integrals and composite external fields. Such factors are incorporated in our framework simultaneously, meaning that the methodology is readily extended to other low-dimensional systems with arbitrary layer number and stacking structure. The predicted results could be validated by magneto-optical measurements[35-42].

Figure 1

Geometric structures of bilayer phosphorene for the (a) top view and (b) side view with various intralayer and interlayer atomic interactions. In the presence of a uniform perpendicular magnetic field, an enlarged unit cell is in rectangular shape as shown in (c). The band structures under various electric fields are presented in (d).

Figure 2

The (a) E-dependent LL spectrum at B = 30 T and the (b–g) corresponding wavefunctions at various E’s.

Figure 3

The x-polarized magneto-optical absorption spectra at B = 30 T for various electric-field strengths of (a) E = 0, (b) 0.1, (c) 0.2, (d) 0.29 and (e) 0.32 (V/Å). The pairs of numbers in (a) indicate the excitation channels.

Figure 7

The E-dependent absorption frequency at B = 30 T under (a) x- and (b) y-polarizations. The color scale represents the intensity of the energy-loss function.

Figure 5

The B-dependent optical absorption frequency at E = 0 under (a) x- and (b) y-polarizations. The threshold excitation channels are indicated by the pairs of numbers. The color scale represents the intensity of the energy-loss function.

Methods

In a single layer of phosphorene, the P atoms lie on the lower (pink circles) and higher (blue circles) planes due to the puckered honeycomb structure, as shown in Fig. 1(a,b). The four atoms in a primitive unit cell (as marked by the dashed green lines in Fig. 1(a)) belong to four different sublattice sites, denoted as A1, A2, B1 and B2. For a N-layer system, there are 4 × N atoms in a unit cell, e.g., 8 atoms for bilayer phosphorene (Fig. 1(b)). The low-lying energy bands are mainly dominated by the atomic interactions of 3 p orbitals. The Hamiltonian of a few-layer phosphorene can be expressed as

Here,  eV is a layer- and sublattice-dependent site energy due to the chemical environment difference between the A and B sublattices in bilayer BP[1]. is the Coulomb potential energy induced by an electric field. They both contribute to the diagonal matrix elements. The summation covers all the lattice sites on two layers. is the annihilation (creation) operator which can destroy (create) an electronic state on the i-th (j-th) lattice site of the l-th (l′-th) layer; i = 1, 2, 3, 4, respectively, correspond to the P atoms on the A1, B1, B2 and A2 sublattices. and represent the intralayer and interlayer hopping integrals, respectively. The effective atomic interactions used in the calculations are denoted by  = 3.665 eV,  eV,  eV,  eV,  eV,  eV,  eV,  eV and  eV[1], as indicated in Fig. 1(a,b).

When bilayer phosphorene is present in a uniform perpendicular magnetic field, the magnetic flux through a unit rectangle is Φ = a1a2B, where a1 = 3.27 Å and a2 = 4.43 Å are lattice constants (the lattice vectors shown by the green arrows in Fig. 1(a)). The vector potential,  = (Bx), can create the Peierls phase of in the hopping integrals, leading to a new period along and thus an enlarged rectangular unit cell with 4R = 4ϕ0/Φ atoms in each phosphorene layer, as illustrated in Fig. 1(c). ϕ0 ( = hc/e = 4.1 × 10−15 T⋅m2) is the magnetic flux quantum; ϕ0/Φ is chosen to be an integer for a model study. The reduced first Brilloun zone has an area of 4π2/a1a2R, in which all the LLs are degenerate for any (k, k) states. For bilayer BP, the magnetic Hamiltonian matrix, which is built from the 8R tight-binding functions, has 8R × 8R elements. This Hamiltonian is a huge Hermitian matrix within achievable experimental field strengths, e.g., the dimension reaches 50400 at B = 10 T.

After the diagonalization of bilayer magnetic Hamiltonian, the LL wave function, with quantum number n, could be expressed aswhere is the tight-binding function localized at the sublattice-dependent lattice sites. / is the amplitude on the -sublattice-dependent lattice site. Specifically, all the amplitudes in an enlarged unit cell could be regarded as the spatial distributions of the sub-envelope functions on the distinct sublattices. Bilayer BP has eight magnetic sub-envelope functions, in which the dominating ones are utilized to characterize the magnetic quantum numbers and the types of LLs. They provide much information for explaining the peculiar LL behaviors, as discussed below. The strong electrically-tunable LL energy spectra exist when B < 60 T. The theoretical framework incorporates the intra-layer and inter-layer atomic interactions, as well as the effects due to geometric structures and external fields, simultaneously. It is suitable for studying the diverse quantization phenomena in layered materials with distinct stacking configurations under uniform/modulated/composite external fields. Moreover, the calculated results are accurate and reliable over a wide energy range.

The electric-dipole optical excitations directly reflect the main characteristics of the electronic properties. The occupied valence LL of n is vertically excited to the unoccupied conduction LL of n in the presence of an electromagnetic wave at zero temperature, being characterized by the excitation channel of Δn = |n−n|. Based on the linear-response Kubo theory, the optical spectral function, without the change of wave vector in the inter-LL transition, is given by

The contributions, which come from the spin, localization-center and k degeneracies, are included in the first factor appearing before the summation. The available transition channels and the excitation intensity depend on the square of the velocity matrix element |M(n, n)|2. The last factor is the joint density-of-states with the broadening parameter γ, which is related to the initial [(n, k)] and final [(n, k)] LLs. Also, is a unit vector of electric polarization and p is momentum. The parallel and perpendicular polarizations, and , are taken into account, respectively. The velocity matrix element M(n, n) is calculated within the gradient approximation and is expressed by the amplitudes ’s(’s) of the tight-binding functions ’s on the distinct sublattices. For , the velocity matrix element is evaluated from

In the gradient approximation, , depending on the intralayer and interlayer hopping integrals combined with the Peierls phases. The quantity M(n, n), which is determined by the relations among eight sub-envelope functions and the sublattice-dependent hopping integrals/Peierls phases, can account for the magneto-optical selection rules. By accurate calculations and detailed examinations of the well-behaved LLs, the product of three terms in Eq. (4) remains the same or changes the sign after the interchange of the initial and final states on the sublattices. Its direct summation leads to a finite (vanishing) velocity matrix element for Δn = odd (even) integers under the x-polarization (discussed latter for Fig. 3) and the opposite behavior is revealed for the y-polarization (Fig. 4).

Figure 4

Similar plot as Fig. 3 for the y-polarization.

Results and Discussion

Landau Level Spectrum

The special lattice structure and complicated hopping integrals generate rich band structures. Bilayer BP has a direct gap of E ~ 1 eV near the Γ point, as illustrated in Fig. 1(d) by the black curves, being smaller than that (~1.6 eV) for the monolayer system. The highly anisotropic energy bands yield the approximately linear and parabolic dispersions along Γ X and Γ Y ( and ), respectively. The energy gap could be reduced considerably by applying an electric field (blue curves). The semiconductor-semimetal transition appears at the strength larger than (V/Å; red curves), for which the valence and conduction bands are transformed into linearly intersecting bands and oscillatory bands along Γ Y and Γ X, respectively. Two split Dirac-cone structures are situated on the right- and left-hand sides of the Γ point (along and ). The extreme points remain at the Γ point, accompanied by two saddle points on the opposite . The electronic states near the Dirac and Γ points will be magnetically quantized into two distinct LL subgroups (discussed in Fig. 2)[31].

All the critical points and constant-energy loops in the energy-wave vector space will dominate the main features of LLs. When E = 0, each LL is four-fold degenerate for each (k, k) state in the presence of spin and localization-center degeneracies. The occupied LLs are asymmetric compared with the unoccupied ones near the Fermi level. The LL energy spectrum is very sensitive to the electric-field strength, as illustrated in Fig. 2(a). The conduction and valence LL energies, respectively, rapidly decline and grow with an increase of E. The lowest unoccupied and the highest occupied LLs, corresponding to band-edge states, merge together at E, in which they are closely related to the magnetic quantization of electronic states near the Dirac points. A larger E can extend the range of linear energy bands and thus double the degeneracy of low-lying LLs (from two neighboring Dirac points in ref.[31]). The LL anti-crossings and crossings occur frequently, since there are two competitive/cooperative LL subgroups initiated from the Dirac and Γ points, respectively[31].

The amplitude, localization center and oscillation mode of the LL wave functions strongly depend on the electric-field strength. For bilayer BP systems, the eight sub-enevlope functions on the distinct sublattices have simple relations. Therefore, is illustrated to understand the critical characteristics and the E-dependence. In the absence of E, and . The quantum number n (n) for each conduction (valence) LL is clearly identified from the number of zero points. For example, the four low-lying conduction/valence LLs have n = 0, 1, 2; 3, as shown in Fig. 2(b). The localization centers are near the 1/2 and 2/2 positions of the enlarged unit cell of the crystal lattice, reflecting the magnetic quantization associated with the Γ point. The former is sufficient for studying the magneto-optical properties. The well-behaved spatial distributions are similar to those in monolayer graphene[43]. In addition, the deeper/higher n LLs might have the side modes of n ± 1, 2; 3 due to the complicated interlayer hopping integrals.

A perpendicular electric field breaks the inversion symmetry, leading to the probability transfer among the distinct sublattices and even the changes of spatial oscillation modes. With a gradual increase of E, the four sub-envelope functions on the first and second layers, respectively, have a simple relation, as revealed in the E = 0 case (Fig. 2(b–d)). Such functions change slowly for the conduction LLs, while the opposite is true for the valence LLs. Their oscillation modes are dramatically changed from n into n + 1 as E gets close to the critical value, e.g., E = 0.29 and 0.3 (Fig. 2(e,f)). And then, the n = 0 and n = 1 LLs are hybridized with each other, so that their wavefunctions present the combination of two oscillation modes (E = 0.32 in Fig. 2(g)). Furthermore, there exist two localization centers on the left- and right-hand sides of the 1/2 position, reflecting the existence of two neighboring Dirac-cone structures[31]. The unusual/perturbed oscillation mode, with a main mode and significant side modes, are also revealed in other low-lying LLs. In addition, the LL splitting will become obvious at deeper/higher energies, in which the Γ-induced LL subgroup comes to exist and exhibits anti-crossing behaviors with the Dirac-point-related one. The various features of sub-envelope functions under distinct electric-field strengths imply the electrically tunable optical selection rules.

Magneto-absorption spectra

The magneto-optical excitations directly reflect the main features of magnetic quantization. The available inter-LL excitation channels, which arise from transitions from the occupied n states to the unoccupied n ones, are denoted as (n, n). They exhibit a lot of delta-function-like absorption peaks, being sensitive to external fields and polarization directions. For the parallel polarization direction, the magneto-absorption symmetric peaks, as shown in Fig. 3(a) at B = 30 T, are dominated by the selection rules of Δn = 1 and 3. The threshold peak comes from the (0, 1) channel. There are many double-peak structures as a result of the asymmetric LL energy spectrum (Fig. 2(a)). Furthermore, the Δn = 1 absorption structures (red circles) might be much higher than the Δn = 3 ones (black circles). Our numerical calculations show that the finite velocity matrix elements related to the former are induced by the symmetric/anti-symmetric superpositions of eight sub-envelopes in the initial and final LL states and the dominating intralayer hopping integral () (Eq. 4). The extra selection rule, such as the latter, is driven by the side modes of the deeper/higher perturbed LLs under the interlayer hopping integrals. Specifically, all the inter-LL absorption peaks present non-uniform intensities, regardless of the selection rules.

The complex relations between the Coulomb potential energies and the intrinsic interactions might dramatically change the magneto-optical selection rules. In addition to the Δn = 1 excitation channels, the Δn = 0 ones (green circles) come to exist quickly as E gradually grows, as clearly demonstrated in Fig. 3(b) for E = 0.1. These two kinds of absorption peaks appear alternatively, in which the neighboring two peaks of the former are merged together and the threshold peak arises from the latter. The Δn = 0 selection rule reflects the fact that an electric field could create an obvious asymmetry between the valence and conduction sub-envelope amplitudes on the same sublattices (Fig. 2(b)), especially for the low-lying LLs. This asymmetry also affects the peak intensities, so that the Δn = 0 peaks are stronger (weaker) than the Δn = 1 ones at lower (higher) frequency. The Δn = 0 channels will become the dominant excitations with the further increase of field strength, e.g., E = 0.2 in Fig. 3(c), a result of the enhanced amplitude asymmetry (Fig. 2(d)). But, as E approaches the critical field, Δn = 1 peaks increase rapidly, compared with the Δn = 0 ones, e.g., E = 0.29 in Fig. 3(d). This is attributed to the E-induced significant side modes in the low-lying valence LLs (Fig. 2(e)). Specially, for E ≥ E, the Dirac-point-created LLs exhibit very strong absorption spectra at lower frequency, e.g., the prominent peaks of Δn = 1, 3; 0 at ω < 0.06 eV under E = 0.32 (the left region of the brown vertical dashed line in Fig. 3(e)). The main reason is that such LLs have double degeneracy and two localization centers with the main and side modes (Fig. 2(g)). Moreover, the Γ-related LLs also make important contributions to the higher-frequency absorption peaks (the right region of the brown vertical dashed line).

The magneto-optical properties are strongly anisotropic even for the degenerate LLs. The spectral intensities in the perpendicular polarization decline significantly as indicated from the comparison between Figs 3(a) and 4(a). The y-polarization velocity matrix element depends on the smaller hopping integrals, but is independent of the largest one (Eq. (4)). This illustrates the drastic changes of peak intensities during the variation of polarization direction. The magneto-optical selection rules become Δn = 2, 0; 4 (blue, pink; yellow circles; discussed in Eq. (4)). The Δn = 1; 3 channels disappear by satisfying a specific relation in the two-sublattice-dependent velocity matrix element. The Δn = 2 and 0 channels, respectively, dominate absorption peaks at lower and higher frequencies. Apparently, the (0, 2) channel creates the threshold peak. When the electric field is below its critical value, only Δn = 2 peaks survive (Fig. 4(b–d)). The Δn = 0 selection rule is greatly suppressed through the strong asymmetry of LL wave-function amplitudes. However, more available excitation channels come to exist for E ≥ E, e.g., the Δn = 1, 3; 0 selection rules at E = 0.32 in Fig. 4(e)). A similar phenomenon, the diversified selection rule, also appears in the x-polarization (Fig. 3(e)), directly reflecting the characteristics of the Dirac-point-induced neighboring LL subgroup and the strong crossings/anti-crossings with the Γ-created LL subgroup.

The dependencies of the magneto-absorption spectra on the magnetic- and electric-field strengths deserve a closer examination. In the absence of E, the peak frequencies, corresponding to the selection rules (Δn = 1; 3 in Fig. 5(a) and Δn = 2, 0; 4 in Fig. 5(b)), rise with increasing B. The B-dependence deviates from the square-root and linear behaviors, being different from those for monolayer graphene and the 2D electron gas. This is closely related to the unusual band structure (Fig. 1(d)). There exist certain important differences between the parallel and perpendicular polarizations in terms of the intensity, frequency and number of absorption peaks. The largest intralayer hooping integral involved in the x-polarization excitations is responsible for the significant anisotropy in the absorption intensity. The threshold frequencies in the x- and y-polarization are, respectively, determined by Δn = 1 and 2, so that optical gaps are lower under the former. The monotonic B-dependence is also revealed in electric fields below the critical one.

On the other side, magneto-optical properties possess the unique B-dependence beyond the critical electric field, as shown in Fig. 6(a,b) at E = 0.32. Two categories of inter-LL excitations, which arise from the Dirac-cone- and Γ-valley-related LL subgroups (Fig. 2(a)), respectively, contribute to lower- and higher-frequency absorption peaks. The LL energies of the former and the latter, respectively, rise and decline with an increase in magnetic field and so do their absorption peak frequencies. The B-dependent absorption frequencies become non-monotonic and abnormal when the anti-crossings of two subgroups frequently appear at sufficiently high magnetic field (B > 30 T). Specifically, the lower absorption frequencies due to the Dirac-cone LLs present the special -dependence for B < 30 T, as observed in monolayer graphene[38,39,44]. The first absorption peak is induced by the occupied/unoccupied LL at the Fermi level and the nearest conduction/valence LL. Its intensity is strongest among all the absorption peaks. The threshold frequency is independent of polarization direction. Furthermore, it increases with B in the monotonic form.

Figure 6

Same plot as Fig. 5, but shown at E = 0.32. The lowest few excitation channels in (a) are marked by the pairs of numbers.

The main features of the magneto-absorption spectra are very sensitive to the electric-field strength. For E < E, the frequencies of absorption peaks decrease rapidly e.g., E < 0.302 in Fig. 7(a,b) at B = 30 T. On the other hand, the unusual magneto-optical properties are revealed in the range of E > E. They present non-monotonic, oscillatory and discontinuous E-dependencies. Since the quantum modes are altered by a sufficiently high electric field, the new/original excitation channels will appear/disappear. Some absorption peaks could only survive for certain ranges of E’s. Obviously, the threshold channel is changed from (0, 0) into (0, 1) [(0, 2) into (0, 1)] near the critical electric field under the x-polarization (the y-polarization). This even leads to the abrupt change of threshold frequency.

The above-mentioned features of the magneto-absorption spectra could be examined by various optical spectroscopies, such as infrared transmission[35-40] and Raman scattering spectroscopies[41,42]. Up to now, experimental measurements have confirmed the rich and diverse magneto-optical properties in graphene-related systems, e.g., Bernal graphite[35-37], AB-stacked few-layer graphene[38,39] and carbon nanotubes[40]. The low-lying absorption peaks in layered graphene have been identified as agreeing with the selection rule of Δn = 1[38,39], in which the square-root and linear dependencies of excitation frequencies on B, respectively, are attributed to the monolayer- and bilayer-like characteristics. The coexistence behavior is also observed in AB-stacked graphite[35-37]. Moreover, the periodic Aharonov-Bohm effect has been verified to exist in cylindrical nanotube systems[40]. The magneto-optical properties in bilayer BP, the selection rule, number, frequency and intensity of absorption peaks, are greatly diversified/enriched by the polarization direction and external fields; that is, they are easily tuned by external factors. The experimental identifications are very useful in understanding the effects due to the geometric symmetry, intrinsic interactions; electric and magnetic fields. It should be noticed that bilayer BP is very different from AB- and AA-stacked graphene[45,46] with respect to the anisotropy/isotropy, selection rules; B- and E-dependencies of magnet-optical spectra. In general, the latter possess almost isotropic excitations, a dominant selection rule with Δn = 1 and a linear or square-root relation between B and frequency of absorption peak.

Conclusions

The generalized tight-binding model, in conjunction with the Kubo formula, has been utilized to investigate the rich and unique magneto-optical properties of bilayer BP. The main features of the LLs, the field-dependent energy spectra and wavefunctions, are well characterized by the oscillation modes of sub-envelope functions on the distinct sublattices. They account for the selection rule, frequency and number of magneto-absorption peaks, strongly depending on the polarization directions; electric and magnetic fields. The predicted results could be verified by magneto-optical spectroscopies, as it has already been done for graphene-related systems[35-42]. The theoretical framework could be further developed to fully understand the diverse magnetic quantization phenomena, especially for the generalization to emergent 2D systems covering few-layer silicene, germanene, tinene, antomonene, bismuthene, etc.

Bilayer BP exhibits the highly anisotropic magneto-absorption spectra, various selection rules and the usual/abnormal field dependencies. The spectral intensity declines obviously during the variation from the x- to y-polarization. This is determined by whether the velocity matrix element is associated with the largest intralayer hopping integral. The selection rules, which come from the well-behaved LLs in the absence of E, are characterized by Δn = odd and even integers for the x- and y-polarizations, respectively. The absorption frequencies deviate from the linear or square-root B-dependence. With the increase of E, the spectral intensities of available channels might be greatly enhanced/suppressed, or the extra selection rules come to exist. The non-monotonic and oscillatory field dependencies are revealed in electric fields beyond the critical one, e.g., the drastic/dramatic changes in threshold frequency/channel near E. Such magneto-optical properties clearly illustrate the close relations among the geometric symmetry, intralayer and interlayer atomic interactions and external fields. There exist important differences between bilayer BP and graphene in the main features of magneto-optical spectra.