Strain-induced gauge and Rashba fields in ferroelectric Rashba lead chalcogenide PbX (X=S, Se, Te) monolayers

Strain-induced gauge and Rashba fields in ferroelectric Rashba lead chalcogenide PbX (X=S, Se, Te) monolayers

Paul Z. Hanakata Department of Physics, Boston University, Boston, MA 02215    A. S. Rodin Centre for Advanced 2D Materials and Graphene Research Centre, National University of Singapore, 6 Science Drive 2, 117546, Singapore    Harold S. Park Department of Mechanical Engineering, Boston University, Boston, MA 02215    David K. Campbell Department of Physics, Boston University, Boston, MA 02215    A. H. Castro Neto Centre for Advanced 2D Materials and Graphene Research Centre, National University of Singapore, 6 Science Drive 2, 117546, Singapore
July 17, 2019

One of the exciting features of two-dimensional (2D) materials is their electronic and optical tunability through strain engineering. Previously we found a new class of 2D ferroelectric Rashba semiconductors PbX (X=S, Se, Te) with tunable spin-orbital properties. In this work, based on our previous tight-binding (TB) results, we derive an effective low-energy Hamiltonian around the symmetry points that captures the effects of strain on the electronic properties of PbX. We find that strains induce gauge fields which shift the Rashba point and modify the Rashba parameter. This effect is equivalent to the application of in-plane magnetic fields. The out-of-plane strain, which is proportional to the electric polarization, is also shown to modify the Rashba parameter. Overall, our theory connects strain and spin-splitting in ferroelectric Rashba materials, which will be important to understand the strain-induced variations in local Rashba parameters that will occur in practical applications.


I Introduction

Monolayers and heterostructures of two-dimensional (2D) materials with spin-orbit interaction offer promise for observing many novel physical effects Geim and Grigorieva (2013); Novoselov et al. (2016); Manchon et al. (2015). In particular, it has been proposed that topological insulators or semiconductors with Rashba interactions coupled with superconductors may host Majorana fermions, which are potential building blocks for topological quantum computers Sau et al. (2010); Fu and Kane (2008).

In addition to 2D materials that exist in the hexagonal phase, such as graphene and the transition metal dichalcogenides (TMDCs), 2D materials with square lattices have been successfully fabricated Chang et al. (2016); Moayed et al. (2017). Recently, the Rashba effect has been observed in thin layers (6–20 nm) of lead sulfide (PbS) Moayed et al. (2017), where an external electric field is used to break the inversion symmetry. However, the spin-splitting is not large. In our previous work based on density functional theory (DFT) calculations, we found that lead chalcogenide monolayers PbX (X=S, Se, Te) have large Rashba coupling  eVÅ in their non-centrosymmetric buckled phase Hanakata et al. (2017). In addition, the spin texture can be switched in a non-volatile way by applying an electric field or mechanical strain, which puts these materials into the family of ferroelectric Rashba semiconductors (FERSCs) Di Sante et al. (2013, 2015). This spin-switching mechanism has recently been observed experimentally in thin films GeTe where the surface is engineered to have either an inward or outward electric polarization Rinaldi et al. (2018).

In reality, monolayers experience strains due to substrates, defects, and so on, where local strains may change the electronic properties of monolayers. Important examples of such effects are pseudo-Landau levels in graphene blisters Levy et al. (2010) and band gap shifts in biaxially strained MoS Castellanos-Gomez et al. (2013). Recently, spatial variations of Rashba coupling due to variations in local electrostatic potentials were reported in InSb Bindel et al. (2016). To date, most theoretical studies of lead chalcogenide monolayers have been based solely on DFT calculations Liu et al. (2015); Wan et al. (2017). However, because DFT is limited to the simulation of small systems, typically several nanometers, it is difficult to model inhomogeneous strains over large spatial areas using DFT.

In this paper, based on our previous tight-binding (TB) model Hanakata et al. (2017); Rodin et al. (2017), we develop a continuum model to predict strain-induced changes in the spin and electronic properties of buckled PbX monolayers. We have also performed DFT calculations to validate our TB predictions. Due to the buckled structure of PbX, the angular dependence becomes important as the relative angle between hybrid orbitals of the top and bottom layer can change substantially Hanakata et al. (2017). We note that some studies on (non-buckled) SnTe and PbX (X=S, Se, Te) rock-salt type materials have incorporated strain effects in the TB, but did not include the changes in hopping parameters due to angle changes Tang and Fu (2014); Barone et al. (2013). In contrast, our TB formulation incorporates the effects due to changes in (i) bond distance and (ii) angle between nearest neighbors as well as (iii) lattice vector deformation.

In the low-energy Hamiltonian, the biaxial (or uniaxial) strains can be described as gauge fields, which are equivalent to, by minimal coupling, the application of in-plane magnetic fields. The out-of-plane strain is directly related to the out-of-plane polarization and this also modifies the Rashba parameter. Within this framework we are able to quantify the Rashba fields in terms of the strain fields.

Figure 1: (a) Schematic top and side views of a buckled monolayer. (b) Undeformed and deformed Brillouin zone as the monolayer is stretched in the and direction. (c) Representative band structures of strained PbS along symmetry points ---- and (d) close to . (e) Relative change in the Rashba parameters obtained from DFT calculations as a function of strain for PbS, PbSe, and PbTe. Energy spin-splitting of PbS for isotropic strains of (f) and (g) . It can be seen that the points are originally located at and shifted closer to the center under a strain of .

Ii Tight-binding

Lead chalcogenide PbX (X=S, Se, Te) consists of two atoms per unit cell, denoted by and atoms, respectively. Lead is a heavy atom (Z(Pb)=82), and it is crucial for creating large spin-orbit interaction (SOI). The schematic top and side views of a buckled lattice are shown in fig. 1(a). is the unit lattice vector and is the vector connecting atom and its neighbor. We denote the relaxed bond length between the neighboring and atoms by , the vector connecting and atoms in the unit cell where , , and is the buckling angle (with corresponding to a flat lattice).

The bands near the Fermi level are mostly composed of and orbitals from both and atoms Hanakata et al. (2017). The bands near the symmetry points can be described within the TB framework including first nearest neighbors and SOI. The full derivation of the TB model can be found in our previous works Hanakata et al. (2017); Rodin et al. (2017), and thus we will only outline the important parts; a more detailed derivation can be found in Appendix C.

For the two atom unit cell shown in Fig. 1(a), the relevant orbital basis involves . To write down the hopping matrix, we use the Slater-Koster matrix elements for the orbitals of neighboring atoms Slater and Koster (1954). As we include the SOI, (where ), we will write our Hamiltonian in angular momentum basis. The dimension of the total Hilbert space is with new basis of , where is the sublattice degree of freedom, is the orbital angular momentum degree of freedom, and is the spin degree of freedom.

We found a Rashba-like dispersion near the and points when the two sublattices are not equivalent Hanakata et al. (2017); Rodin et al. (2017). In this paper, we develop a continuum strain model describing changes in the Rashba dispersion near the point, and thus the Hamiltonian is expanded around the point . Exactly at [], the Hamiltonian decomposes into several uncoupled blocks and the wave function of the conduction band is given by , with , , and being real numbers Hanakata et al. (2017); Rodin et al. (2017). The Hamiltonian for the valance band can be obtained by interchanging and .

Projecting the Hamiltonian onto the conduction band subspace we obtain the effective Rashba-like Hamiltonian


where is the momenta, , is the Rashba parameter, and . The coefficients can be obtained from the DFT calculations. Since we know the buckling angle we can can evaluate . All of the relevant (unstrained) parameters are tabulated in Appendix A.1.

Iii Strain-induced gauge fields

Figure 2: Schematic changes in the Rashba dispersions due to (a) in-plane strains and (b) out-of-plane strains. The linear Rashba dispersions at the for unstrained systems are colored blue. Under positive in-plane strains, the Rashba points shift closer to and the strength of Rashba parameters decrease (smaller slope) with increasing strains. On the other hand, under out-of-plane strain, the strength of Rashba parameters increases with increasing uniaxial out-of-plane strain while the Rashba points do not shift.

Since the SOI is independent of lattice distortions, in this derivation we will focus on the spinless Hamiltonian and then reintroduce the spin terms. We will focus on the conduction band only, as the changes in valence band should be similar.

Under deformation a vector connecting two points in a unit cell can be approximated as , where is the displacement vector, and . In this work we focus on deformation that does not involve local rotation . Similarly, between two lattice vectors .

Alterations in bond distance will result in changes in the hopping energies. Since studies of lead chalcogenides under strain are very limited, we follow the Wills-Harrison’s argument Harrison (2004) and assume that the hopping energy . Similar considerations also have been used for strained TMDCs Cazalilla et al. (2014); Rostami et al. (2015); Fang et al. (2017) and phosphorene Jiang and Park (2015); Taghizadeh Sisakht et al. (2016). Note that the hopping matrix derived from Slater-Koster has angular dependence and these relative angles should change due to strain. Assuming the hopping matrix depends on bond distance only, the modified hopping parameter, in terms of the strain tensor , is  Cazalilla et al. (2014); Rostami et al. (2015). This approximation is also the case for graphene, where the hopping modulation is approximated as . In particular, this approximation works well for flat graphene under strain because the angle between orbitals does not change. The angular dependence becomes more important when deformations, such as nanobubbles and kirigami patterns, create large curvature (bending) Qi et al. (2014a, b). In buckled lead chalcogenides, however, the relevant hopping terms for the Rashba dispersion depend on the buckling angle even in the simple case of biaxial strains Hanakata et al. (2017). Thus we will include this angular dependence, and we will show that this is important to capture the changes in Rashba coupling with uniaxial strain.

Let the unstrained vector connecting an atom and its neighbor be defined as and the equilibrium distance . Here we show the derivation for , while the others can be found by following the same procedure. We assume and we expect  Harrison (2004). In Cartesian coordinates the strained hopping is given by , and by Taylor expansion we obtain,


Within the strain approximation . If we alter only the bond distance while keeping the angle constant, we will get the same expression as above when angular effects are assumed to be negligible.

The interlattice-spinless Hamiltonian in reciprocal space can be written as


where is the sum over nearest neighbor pairs and . The first term is the unstrained Hamiltonian, is the correction due to lattice deformation, and is the correction from the altered hopping parameter due to changes in both the interatomic distance and angle between orbitals.

Iv Homogenous isotropic strains

We start with a simple deformation with no shear We will focus on the matrix elements that are relevant to the conductions band, such as and . In the angular momentum basis, the correction from and at is given by


where , and . Note that are the unstrained geometrical and hopping parameters. is independent of the direction strains (e.g ) because the lattice vector and are two-dimensional. Because of the symmetry of , we found that the first correction at due to bond alterations is first order in and momentum . In graphene, the first correction from hopping modulation that is linear in (but not proportional to ) is not zero Guinea et al. (2010); Mañes et al. (2013); Masir et al. (2013). We have to include the contributions of up to first order in as well because in (-dependent term) we keep terms up to first order in and .

To obtain we will consider an isotropic strain . Notice that the change in low-energy Hamiltonian of Eq. 1 due to and at can be written as gauge potentials,


where where we have used to simplify and is the unstrained Rashba parameter.

and the second term of are proportional to . This modifies the strength of Rashba parameter . This alteration in the Rashba term is similar to the modification of Fermi velocity in graphene de Juan et al. (2012); Mañes et al. (2013); Masir et al. (2013).

We next present our DFT results to validate our TB predictions. Details of DFT calculations and the unstrained geometrical parameters of buckled PbS, PbSe, and PbTe can be found in Appendix A.1. Strains are applied to the relaxed buckled phase. In order to find the effects that come from changes in bond distance only, we deformed the monolayer in the DFT simulations by changing the bond distance while keeping the angle constant. The lattice vectors and atomic positions are not relaxed under this deformation. The Rashba parameters are obtained by taking the derivative of the energy dispersion in the vicinity of the point, . Under isotropic deformations, we found that at decreases with increasing strain (weakening of the hopping interaction), as expected from Eq. 5, shown in fig. 1(c)-(e). A direct comparison between DFT results and TB with strain-included allows us to extract . By fitting DFT data points to a straight line, we obtained for PbS, PbSe, and PbTe, respectively (fig. 1(e)). We see that the value of would be different if the lattice deformation correction was not included.

As we stretch the lattice, the Brillouin zone (BZ) will shrink, and the corner of the BZ ( point) will shift as , where is the undeformed lattice constant. For positive strains, the point shifts towards the point (relative to the undeformed BZ), shown in fig. 1(b). In our modified TB model, the point is displaced due to the first term of the lattice deformation correction (see Eq. 5). The momentum shifts due to lattice deformations are also found in graphene Kitt et al. (2012). The changes in Rashba dispersion and its locations due to strains are illustrated in fig. 2.

To show the momentum shifts relative to the undeformed (reference) state, we plot the energy spin-splitting at the conduction band of PbS obtained from the DFT results as a function of , shown in fig. 1(f) and (g). Note that momenta are in units of . Originally the points are located at and are shifted closer to () when an isotropic strain of is applied. The momentum shift is linear with strains , consistent with several previous works Kitt et al. (2012); Masir et al. (2013). This Rashba-point shift due to strains is equivalent to applying in-plane-magnetic fields to the system,


where , , and is the Bohr magneton. For completeness the derivation of Eq. 6 is included in Appendix D. As an illustration, we can choose an external field of , upon which the in-plane magnetic field is given by . Since the Bohr magneton is small, in order to get a similar effect of 2% strain using magnetic fields, one has to apply external magnetic fields with an approximate strength of Tesla (by Eq. 5 and Eq 6).

Figure 3: (a) Out-of-plane polarization as a function of out-of-plane strain . (b) Linear relationship between and which is consistent with TB predictions. (c) Rashba parameter as a function of . All data points are obtained from the DFT calculations.

V Electric polarization and Rashba field

Proposals have been made to change the spin texture (i.e. sign of ) by changing the electric polarization Di Sante et al. (2013); Liebmann et al. (2016); Leppert et al. (2016); Liu et al. (2013). Rinaldi et al. found that the spin-texture in FERSC GeTe films indeed depends on the locations of the atoms on the surface, which dictate the direction of the electric polarization Rinaldi et al. (2018). In DFT simulations of SnTe thin films, which have a structure similar to PbX, it also has been shown that near the vacuum one of the atomic species buckles outward while the other species buckles inward Qian et al. (2015). While the proportionality between Rashba parameter and spontaneous electric polarization is well known, it will be useful to understand this mechanism in PbX from a microscopic view, where the changes in Rashba parameters can be understood in terms of interactions between atoms and the external applied strains. We will show that our strain-dependent TB model captures how the out-of-plane strain, which is proportional to the out-of-plane polarization, modifies the Rashba fields.

By the modern theory of polarization, the electric polarization is given byKing-Smith and Vanderbilt (1993) , where is the ionic charge plus the core electrons, is the position of ions, is the unit cell volume, is the elementary charge, is the valence band index, is the wave vector, and is the electronic wave function. The first term is the contribution from core electrons and ions, and the second term is the electronic contribution defined as the adiabatic flow of current, which can be calculated from the Berry phase (BP) King-Smith and Vanderbilt (1993). The spontaneous polarization is calculated by taking the difference between the polarization of the polar (buckled) state and the non-polar (reference) state, . We estimate the thickness to be 0.5 nm in order to compare the polarizations to typical bulk ferroelectrics. Details can be found in Appendix B. In the DFT simulations we distort the ions in the direction while keeping the in-plane lattice vectors fixed at the relaxed buckled values. We report only the spontaneous polarizations of PbS and PbSe, as PbTe is metallic Hanakata et al. (2017). A modified Berry phase calculation is needed to evaluate polarization of ferroelectric metals Filippetti et al. (2016); however this is beyond the scope of our present study.

From the DFT results we found that the core electronic plus ionic and the electronic contribution (BP) are proportional to the distance between Pb and X (X=S, Se) in the direction (plotted in Appendix B). This gives a proportionality between and , as shown in fig. 3(a). Compressing the monolayer in the with strain results in a decrease in , shown in fig. 3(b). This is opposite to the case of isotropic deformation (see fig. 1(e)). This result is consistent with TB predictions. In the previous discussion, we found that increasing bond distance ( generally weakens the hopping interaction and thus decreases . Using relaxed geometrical parameters (i.e buckling angle ) and from Eq. 4, is expected to decrease with compressive strain in the as is negative. We also want to note that there is no gauge-field since is two-dimensional, and thus is not shifted. The changes in Rashba dispersion and its locations due to out-of-plane strains are illustrated in fig. 2(b). Notice that not including the angular dependence in the hopping correction will not capture this effect. The inclusion of the angular dependence is particularly important for the PbX monolayer due to its buckled nature. Overall, this suggests that the out-of-plane internal electric polarization acts as an in-plane gauge field in the low-energy Hamiltonian. Assuming small strains, we found that . This result is important as it establishes a direct relationship between the Rashba field and the out-of-plane polarization which is also proportional to the out-of-plane strain . Recently, several works have also studied strain-induced piezoelectricity in boron nitride Droth et al. (2016) and TMDCs Rostami et al. (2018). Several experimental works use out-of-plane magnetic fields (parallel to the polar axis of Rashba materials) to measure the Rashba parameter as the Landau level spectrum changes with the strength of the Rashba parameter Bindel et al. (2016); Bordács et al. (2013). One could also use this experimental approach to detect variations in the Rashba parameter in PbX due to out-of-plane strains.

Vi Conclusions

We have developed a TB model where the electronic changes in PbX can be described within continuum mechanics. We found the scaling exponent that modifies the hopping parameter to be . In the low-energy Hamiltonian, the effect of strains can be described as gauge fields, which are equivalent to, by minimal coupling, application of an in-plane magnetic field. Our theory describes how the location of the Rashba point and the strength of the Rashba field can be engineered by applying strains. The out-of-plane strain in particular is directly related to the out-of-plane polarization. Within this framework we are able to understand the connection between the Rashba and ferroelectricity.

Our strain-dependent TB model should be applicable for calculating the effects of inhomogeneous strain on the spatially-resolved Rashba fields over a large region, whereas this calculation would not be feasible within a reasonable time using a DFT approach. Employing classical atomistic simulations (e.g. molecular dynamics) together with strain-dependent TB will be an efficient tool for studying larger and more realistic systems with strain modulation due to substrates, indentors Levy et al. (2010); Castellanos-Gomez et al. (2013); Georgi et al. (2017) or geometrical cuts Qi et al. (2014b); Hanakata et al. (2016a).This will open possibilities of using lead chalcogenides for strain and electric-controlled spintronic devices.


  • Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  • Novoselov et al. (2016) K. Novoselov, A. Mishchenko, A. Carvalho,  and A. C. Neto, Science 353, aac9439 (2016).
  • Manchon et al. (2015) A. Manchon, H. C. Koo, J. Nitta, S. Frolov,  and R. Duine, Nature materials 14, 871 (2015).
  • Sau et al. (2010) J. D. Sau, R. M. Lutchyn, S. Tewari,  and S. Das Sarma, Phys. Rev. Lett. 104, 040502 (2010).
  • Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  • Chang et al. (2016) K. Chang, J. Liu, H. Lin, N. Wang, K. Zhao, A. Zhang, F. Jin, Y. Zhong, X. Hu, W. Duan, Q. Zhang, L. Fu, Q.-K. Xue, X. Chen,  and S.-H. Ji, Science 353, 274 (2016).
  • Moayed et al. (2017) M. M. R. Moayed, T. Bielewicz, M. S. Zöllner, C. Herrmann,  and C. Klinke, Nature Communications 8, 15721 (2017).
  • Hanakata et al. (2017) P. Z. Hanakata, A. S. Rodin, A. Carvalho, H. S. Park, D. K. Campbell,  and A. H. Castro Neto, Phys. Rev. B 96, 161401 (2017).
  • Di Sante et al. (2013) D. Di Sante, P. Barone, R. Bertacco,  and S. Picozzi, Advanced Materials 25, 509 (2013).
  • Di Sante et al. (2015) D. Di Sante, A. Stroppa, P. Barone, M.-H. Whangbo,  and S. Picozzi, Physical Review B 91, 161401 (2015).
  • Rinaldi et al. (2018) C. Rinaldi, S. Varotto, M. Asa, J. Sławińska, J. Fujii, G. Vinai, S. Cecchi, D. Di Sante, R. Calarco, I. Vobornik, G. Panaccione, S. Picozzi,  and R. Bertacco, Nano Letters 18, 2751 (2018), pMID: 29380606, .
  • Levy et al. (2010) N. Levy, S. Burke, K. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. C. Neto,  and M. Crommie, Science 329, 544 (2010).
  • Castellanos-Gomez et al. (2013) A. Castellanos-Gomez, R. Roldán, E. Cappelluti, M. Buscema, F. Guinea, H. S. van der Zant,  and G. A. Steele, Nano letters 13, 5361 (2013).
  • Bindel et al. (2016) J. R. Bindel, M. Pezzotta, J. Ulrich, M. Liebmann, E. Y. Sherman,  and M. Morgenstern, Nature Physics 12, 920 (2016).
  • Liu et al. (2015) J. Liu, X. Qian,  and L. Fu, Nano letters 15, 2657 (2015).
  • Wan et al. (2017) W. Wan, Y. Yao, L. Sun, C.-C. Liu,  and F. Zhang, Advanced Materials 29, 1604788 (2017), 1604788.
  • Rodin et al. (2017) A. S. Rodin, P. Z. Hanakata, A. Carvalho, H. S. Park, D. K. Campbell,  and A. H. Castro Neto, Phys. Rev. B 96, 115450 (2017).
  • Tang and Fu (2014) E. Tang and L. Fu, Nature Physics 10, 964 (2014).
  • Barone et al. (2013) P. Barone, D. Di Sante,  and S. Picozzi, Physica status solidi Rapid Research Letters 7, 1102 (2013).
  • Slater and Koster (1954) J. C. Slater and G. F. Koster, Physical Review 94, 1498 (1954).
  • Harrison (2004) W. A. Harrison, Elementary Electronic Structure: Revised (World Scientific Publishing Company, 2004).
  • Cazalilla et al. (2014) M. A. Cazalilla, H. Ochoa,  and F. Guinea, Phys. Rev. Lett. 113, 077201 (2014).
  • Rostami et al. (2015) H. Rostami, R. Roldán, E. Cappelluti, R. Asgari,  and F. Guinea, Physical Review B 92, 195402 (2015).
  • Fang et al. (2017) S. Fang, S. Carr, J. Shen, M. A. Cazalilla,  and E. Kaxiras, arXiv preprint arXiv:1709.07510  (2017).
  • Jiang and Park (2015) J.-W. Jiang and H. S. Park, Phys. Rev. B 91, 235118 (2015).
  • Taghizadeh Sisakht et al. (2016) E. Taghizadeh Sisakht, F. Fazileh, M. H. Zare, M. Zarenia,  and F. M. Peeters, Phys. Rev. B 94, 085417 (2016).
  • Qi et al. (2014a) Z. Qi, A. L. Kitt, H. S. Park, V. M. Pereira, D. K. Campbell,  and A. H. Castro Neto, Phys. Rev. B 90, 125419 (2014a).
  • Qi et al. (2014b) Z. Qi, D. K. Campbell,  and H. S. Park, Phys. Rev. B 90, 245437 (2014b).
  • Guinea et al. (2010) F. Guinea, M. Katsnelson,  and A. Geim, Nature Physics 6, 30 (2010).
  • Mañes et al. (2013) J. L. Mañes, F. de Juan, M. Sturla,  and M. A. H. Vozmediano, Phys. Rev. B 88, 155405 (2013).
  • Masir et al. (2013) M. R. Masir, D. Moldovan,  and F. Peeters, Solid State Communications 175, 76 (2013).
  • de Juan et al. (2012) F. de Juan, M. Sturla,  and M. A. H. Vozmediano, Phys. Rev. Lett. 108, 227205 (2012).
  • Kitt et al. (2012) A. L. Kitt, V. M. Pereira, A. K. Swan,  and B. B. Goldberg, Phys. Rev. B 85, 115432 (2012).
  • Liebmann et al. (2016) M. Liebmann, C. Rinaldi, D. Di Sante, J. Kellner, C. Pauly, R. N. Wang, J. E. Boschker, A. Giussani, S. Bertoli, M. Cantoni, L. Baldrati, M. Asa, I. Vobornik, G. Panaccione, D. Marchenko, J. Sánchez-Barriga, O. Rader, R. Calarco, S. Picozzi, R. Bertacco,  and M. Morgenstern, Advanced Materials 28, 560 (2016).
  • Leppert et al. (2016) L. Leppert, S. E. Reyes-Lillo,  and J. B. Neaton, The journal of physical chemistry letters 7, 3683 (2016).
  • Liu et al. (2013) Q. Liu, Y. Guo,  and A. J. Freeman, Nano Letters 13, 5264 (2013).
  • Qian et al. (2015) X. Qian, L. Fu,  and J. Li, Nano Research 8, 967 (2015).
  • King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  • Filippetti et al. (2016) A. Filippetti, V. Fiorentini, F. Ricci, P. Delugas,  and J. Íñiguez, Nature communications 7, 11211 (2016).
  • Droth et al. (2016) M. Droth, G. Burkard,  and V. M. Pereira, Phys. Rev. B 94, 075404 (2016).
  • Rostami et al. (2018) H. Rostami, F. Guinea, M. Polini,  and R. Roldán, npj 2D Materials and Applications 2, 15 (2018).
  • Bordács et al. (2013) S. Bordács, J. G. Checkelsky, H. Murakawa, H. Y. Hwang,  and Y. Tokura, Phys. Rev. Lett. 111, 166403 (2013).
  • Georgi et al. (2017) A. Georgi, P. Nemes-Incze, R. Carrillo-Bastos, D. Faria, S. Viola Kusminskiy, D. Zhai, M. Schneider, D. Subramaniam, T. Mashoff, N. M. Freitag, M. Liebmann, M. Pratzer, L. Wirtz, C. R. Woods, R. V. Gorbachev, Y. Cao, K. S. Novoselov, N. Sandler,  and M. Morgenstern, Nano Letters 17, 2240 (2017), pMID: 28211276, .
  • Hanakata et al. (2016a) P. Z. Hanakata, Z. Qi, D. K. Campbell,  and H. S. Park, Nanoscale 8, 458 (2016a).
  • Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari,  and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009).
  • Perdew et al. (1996) J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  • Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  • Kobayashi (2015) K. Kobayashi, Surface Science 639, 54 (2015).
  • Hanakata et al. (2016b) P. Z. Hanakata, A. Carvalho, D. K. Campbell,  and H. S. Park, Physical Review B 94, 035304 (2016b).
  • Fei et al. (2016) R. Fei, W. Kang,  and L. Yang, Phys. Rev. Lett. 117, 097601 (2016).
  • Wang and Qian (2017) H. Wang and X. Qian, 2D Materials 4, 015042 (2017).
  • Stroppa et al. (2015) A. Stroppa, C. Quarti, F. De Angelis,  and S. Picozzi, The journal of physical chemistry letters 6, 2223 (2015).
P.Z.H. developed the theory, wrote the paper, and performed the DFT calculations. A.S.R. contributed to the analytical work. H.S.P., D.K.C., and A.C.H.N. supervised the research. We thank Vitor M. Pereira for useful comments and discussions. P.Z.H., H.S.P., and D.K.C. acknowledge support by the Boston University Materials Science and Engineering Innovation Grants. P.Z.H, H.S.P., and D.K.C are grateful for computing resources provided by theBostonUniversity Shared Computing Cluster.

Appendix A Methods

a.1 Computational Details

To validate our tight-binding predictions we performed density functional theory (DFT) calculations implemented in the Quantum ESPRESSO package Giannozzi et al. (2009). We employed projector augmented-wave (PAW) type pseudopotentials with Perdew-Burke-Ernzerhof (PBE) within the generalized gradient approximation (GGA) for the exchange and correlation functional with Perdew et al. (1996). The Kohn-Sham orbitals were expanded in a plane-wave basis with a cutoff energy of 100 Ry and a charge density cutoff of 200 Ry. The cutoff was chosen following the suggested minimum cutoff in the pseudopotental file. A -point grid sampling was generated using the Monkhorst-Pack scheme with 16161 points Monkhorst and Pack (1976). A vacuum of 20 Å was used. The relaxed structures of PbS, PbSe, and PbTe were obtained by relaxing the ionic positions and the lattice vectors. A convergence threshold on total energy of eV and a convergence threshold on forces of 0.005 eV/ were chosen. Lattice vectors are relaxed until the stress is less than 0.01 GPa. Our first-principles calculations show that the buckled phase of the PbX monolayer is more stable than the centrosymmetric planar phase Hanakata et al. (2017), consistent with other DFT studies Wan et al. (2017); Kobayashi (2015). Detailed discussions on the bistable nature, ferroelectric properties and orbital-spin texture properties of lead chalcogenides can be found in our previous paper Hanakata et al. (2017). In the current work, the deformations (atomic distortions) are applied to the optimized buckled structure.

We used a finer grid for band structure calculations with the spin-orbit interaction included. We have tried several large numbers of points and found that a grid of 100 points between two symmetry points (e.g between and ) is enough to obtain the Rashba parameter at the point Hanakata et al. (2017). A regular grid of 40401 was used for the surface plot of the energy spin splitting.

Here we tabulate the optimized (relaxed) geometrical parameters of buckled PbX (X=S, Se, and Te) monolayers in table 1. The Rashba parameters are obtained by taking the derivative of energy dispersion near the point. The orbital coefficients are obtained by projecting the wave functions into the atomic orbital basis. The unstrained values of and are tabulated in table 2. From the table it can be seen that the wave functions are mostly composed of in-plane and out-of-plane of orbitals of Pb and an in-plane orbital of the chalcogen X (X=S, Se, Te).

(Å) (Å) (Å)
PbS 3.74 21.6 1.04 2.84
PbSe 3.82 24.3 1.22 2.96
PbTe 4.01 26.3 1.40 3.16
Table 1: Relaxed lattice constant , buckling angle , buckling height , nearest-neighbor bond distance .
(eVÅ) (eV)
PbS 0.305 0.534 0.115 3.40 5.36
PbSe 0.272 0.549 0.137 3.37 4.28
PbTe 0.286 0.522 0.130 3.18 3.83
Table 2: Rashba parameters , projected wave functions coefficients , , obtained from DFT and .

Appendix B Electric polarization

We used the modern theory of polarization King-Smith and Vanderbilt (1993) to calculate the spontaneous polarization implemented in the Quantum ESPRESSO package Giannozzi et al. (2009). The electric polarization is calculated via Berry phase calculation King-Smith and Vanderbilt (1993), which is given by


where is the ionic charge plus the core electrons, is the position of ions, is the unit cell volume, is the elementary charge, is the valence band index, is the wave vector, and is the electronic wave function. The first term is the contribution from core electrons and ions, and the second term is the electronic contribution defined as the adiabatic flow of current which can be calculated from the Berry connection King-Smith and Vanderbilt (1993).

The spontaneous polarization is calculated by taking the difference between the polarization of the polar (buckled) state and the non-polar (reference) state, . To find the polarization at different heights, we change the out-of-plane distance between the Pb and X (X=S, Se) atom while keeping the in-plane lattice vectors fixed at the optimized buckled values. It is a common practice to use a value on the order of the bulk lattice constant (0.5–1 nm) to estimate the monolayer thickness in order to compare the polarizations of monolayers to the typical bulk ferroelectrics Hanakata et al. (2016b); Fei et al. (2016); Wang and Qian (2017). In this current work we estimate the thickness to be 0.5 nm. In Quantum ESPRESSO, spontaneous polarization with spin-orbit included can be calculated using norm conserving pseudopotentials. A difference of 0.03 is found when spin-orbit interaction is included. Thus, to save computational time we only report spontaneous polarization without inclusion of the spin-orbit interaction. This small difference has also been reported previously Stroppa et al. (2015); Leppert et al. (2016). In figure. 4 we plot the polarization from the ionic plus core electron contribution, and the electronic contribution, from the Berry phase calculation, scaled by their values at zero strain as a function of distance between the Pb and S atom in the direction.

Figure 4: From the DFT results we found that the ionic plus core electronic and the electronic (by Berry phase calculation) contributions are proportional to the distance between Pb and X (X=S, Se) in the direction.

Appendix C Tight binding

The lead chalcogenide monolayer has two atoms per unit cell (). Based on density functional theories, the relevant orbitals near the valence and conduction bands are and orbitals. The wave function of sublattice then can be written as


where is the lattice vector, is a wave vector, is the basis wave function . Including only nearest neighbor hopping the spinless Hamiltonian can be written as


where runs over the onsite cell and the nearest neighboring cells. creates an electron in the unit cell with atomic orbital . We can write this more compactly as


where (the onsite term) is given by


To write down the hopping matrix, we use the following Slater-Koster matrix elements for the orbitals of neighboring atoms Slater and Koster (1954):


Here, is the orientation of the th orbital and is the unit vector pointing from atom 1 to atom 2. If we include up to first nearest neighbors only we can write the inter-lattice hopping matrix as


where . The momentum and . To keep the expression more compact, we have introduced . In addition, since the and species are not necessarily the same, we have two quantities of the form.

While it is convenient to use and orbitals to write down the hopping matrix, since we are interested in including SOI in our model, it is helpful to go to a basis which is more natural for the angular momentum operators:


where the first number represents the orbital momentum quantum number and the second one is its projection along the direction. This basis change does not alter the and matrices. The inter-lattice hopping portion of the Hamiltonian, on the other hand, becomes


From here we write , where is a matrix projector from the orbital basis to the angular momentum basis.

To include the SOI, we use the standard form describing the spin-orbit coupling arising from the interaction with the nucleus:


where is either Pb or X (X=S, Se, Te). The last term modifies the diagonal elements of the self-energy for by adding (subtracting) if and point in the same (opposite) direction. The first tem couples with and with with the coupling strength .

The total Hamiltonian can then be written as


c.1 Point

We first look around the point . To the leading order in , the hopping matrix is given by,


where is the angle measured from the direction. At (), the Hamiltonian decomposes into several uncoupled blocks with the corresponding bases:


where labels the sublattices and the middle ket denotes the spin state. Using the direct sum notation, we can write down the total Hamiltonian as .

From , we see that for a given , the eigenstates are spin-degenerate. The degeneracy becomes four-fold if the atoms of sublattices and are the same, leading to . Equation (18) shows that at finite there is no coupling between the degenerate states that is linear in momentum. This means that the bands composed of orbitals have local extrema at the point.

Next, we turn to from Eq. (19). Just like for , the bands are doubly or four-fold degenerate depending on whether the sublattices are composed of the same atomic species. Without making assumptions about the lattice composition, the general form of the degenerate states is


with , , and real. At finite ,


where is the two-dimensional Levi-Civita symbol. This coupling between the degenerate states leads to an effective Rashba-like Hamiltonian:


or in the matrix form


We use values of , , and obtained from DFT results. To give better physical pictures of these coefficients, we will solve the Hamiltonian Eq. 19. We treat the spin orbit interaction (SOI) as perturbations and we will assume that where is the index denoting Pb with strong SOI and denotes weak SOI of chalcogen atom. Focusing on , Eq. 19 becomes


and the perturbation


We first solved Eq. 24 to find the eigenvalues and eigenvectors and used first order perturbation theory to obtain the corrections to the eigenvectors. Using MATHEMATICA, we found to the first order in that


Recall that we defined Rashba parameter . From Eq. 26 we see that weakly depends on strains. For this reason, in the main text we assumed and are constant and the corrections to come mostly from and .

Appendix D Magnetic Field

Let us now try to include external fields to the system. The magnetic field can be included via the Peierls substitution so that , where is the vector potential. In addition, applying an external magnetic field leads to the interaction of the electron angular momentum with the field.

The total magnetic moment of an electron is given by


so that


Setting gives