⇦Chinese Physics Letters, 2024, Vol. 41, No. 3, Article code 037102 Anisotropic Band Evolution of Bulk Black Phosphorus Induced by Uniaxial Tensile Strain Yafeng Deng (邓亚丰)1, Yilin Zhang (张艺琳)2, Yafei Zhao (赵亚飞)3*, Yongkang Xu (徐永康)1, Xingze Dai (代兴泽)1, Shuanghai Wang (王双海)1, Xianyang Lu (陆显扬)1, Yao Li (黎遥)1, Yongbing Xu (徐永兵)1*, and Liang He (何亮)1* Affiliations 1National Key Laboratory of Spin Chip and Technology, School of Electronic Science and Engineering, Nanjing University, Nanjing 210023, China 2College of Optical and Electronic Technology, China Jiliang University, Hangzhou 310018, China 3School of Physics and Engineering, Henan University of Science and Technology, Luoyang 471023, China Received 21 December 2023; accepted manuscript online 1 February 2024; published online 12 March 2024 *Corresponding authors. Email: zhaoyafei@haust.edu.cn; ybxu@nju.edu.cn; heliang@nju.edu.cn Citation Text: Deng Y F, Zhang Y L, Zhao Y F et al. 2024 Chin. Phys. Lett. 41 037102    Abstract We investigate the anisotropic band structure and its evolution under tensile strains along different crystallographic directions in bulk black phosphorus (BP) using angle-resolved photoemission spectroscopy and density functional theory. The results show that there are band crossings in the Z–L (armchair) direction, but not in the Z–A (zigzag) direction. The corresponding dispersion-$k$ distributions near the valence band maximum (VBM) exhibit quasi-linear or quadratic relationships, respectively. Along the armchair direction, the tensile strain expands the interlayer spacing and shifts the VBM to deeper levels with a slope of $-16.2$ meV/% strain. Conversely, the tensile strain along the zigzag direction compresses the interlayer spacing and causes the VBM to shift towards shallower levels with a slope of 13.1 meV/% strain. This work demonstrates an effective method for band engineering of bulk BP by uniaxial tensile strain, elucidates the mechanism behind it, and paves the way for strain-regulated optoelectronic devices based on bulk BP.
cpl-41-3-037102-fig1.png
cpl-41-3-037102-fig2.png
cpl-41-3-037102-fig3.png
cpl-41-3-037102-fig4.png
cpl-41-3-037102-fig5.png
 DOI:10.1088/0256-307X/41/3/037102 © 2024 Chinese Physics Society Article Text Black phosphorus is a layered van der Waals crystal, shares structural similarities with graphite, and belongs to the orthorhombic crystal system. The exceptional electrical properties of black phosphorus (BP) are exemplified by its significantly high electron mobility[1,2] and high on/off ratio of up to $10^{5}$ in transistor applications.[3-6] In the field of optoelectronics, the direct bandgap of BP can be regulated over a wide range (0.3–2.0 eV) via band engineering,[7-13] and results in effective coupling with various optical wavelengths. The layered BP also demonstrates impressive mechanical resilience, with a critical tensile strain of up to 32% in theory.[14,15] This provides the foundation for research into strain-induced band engineering. Furthermore, due to the in-plane anisotropy of its crystal structure, BP exhibits substantial directional dependence in electrical and optoelectronic properties.[16-18] Numerous reports have explored the strain regulation of the band structure in BP.[9,14,19,20] However, the majority of these studies have primarily concentrated on ultra-thin layers,[21-23] and have used optical analysis to indirectly glean information about the bandgap.[24-27] Thus, there remains a dearth of direct observation of the strain-induced band regulation, especially in bulk BP. In this study, strain-ARPES (angle-resolved photoemission spectroscopy) is employed to study the band structure of bulk BP along the armchair and zigzag directions, respectively. For the first time, a vivid and comprehensive illustration of the anisotropic band structure and its evolution along different crystallographic directions in bulk BP are obtained. More importantly, upon the uniaxial tensile strain along the armchair or the zigzag directions, the valence band maximum (VBM) shifts to deeper or shallower levels, resulting in a VBM shift rate of $-16.2$ meV/% strain or 13.1 meV/% strain, respectively. This anisotropic response is attributed to expansion or compression of the interlayer spacing under different strains. As illustrated in Fig. 1(a), bulk BP consists of numerous wrinkled single layers, held together by van der Waals interactions between the layers. Within each layer, a wrinkled honeycomb-like structure is bonded through $sp^{2}$ and $sp^{3}$ hybridized P atoms. As a two-dimensional layered semiconductor material, BP exhibits pronounced intralayer anisotropy along the armchair [Fig. 1(b)] and zigzag [Fig. 1(c)] directions, making its band structure more easily tunable compared to the in-plane isotropic graphene.[28,29] Simultaneously, in the crystal lattice of bulk BP, the armchair and the zigzag directions correspond to the Z–L and the Z–A directions, respectively, within the momentum space.[25,30,31] This convenient correspondence facilitates our understanding of the relationship between the in-plane lattice structure and the band characteristics. Our initial investigation focused on the band dispersion of bulk BP along both the Z–L–Z (armchair direction) and Z–A–Z (zigzag direction), as shown in Fig. 1(d). Figure 1(e) presents the results of ARPES experiments of bulk BP; the distinct and complete band structure reveals pronounced anisotropic band dispersion along the two directions. Remarkably, the band dispersion along the Z–L direction exhibits a distinctive “bowknot-like” cross-band structure, while the band dispersion along the Z–A direction tends to assume a flatter and non-crossing band distribution.
cpl-41-3-037102-fig1.png
Fig. 1. (a) The crystal structure of bulk BP. (b) Side view. (c) Top view. The lattice constants $a$, $b$, and $c$ follow the zigzag, out of plane, and armchair directions, respectively. (d) The first Brillouin zone of bulk BP. The direct bandgap is located at Z. The anisotropic band structure of bulk BP: (e) experiment, (f) calculation. The band dispersion along the Z–L direction exhibits a distinctive “bowknot-like” cross-band structure, while the band dispersion along the Z–A direction tends to assume a flatter and non-crossing band distribution.
cpl-41-3-037102-fig2.png
Fig. 2. (a) The phonon vibration modes in BP crystal. The Raman shift of bulk BP regulated by uniaxial tensile strain along [(b), (c)] the armchair direction and [(d), (e)] the zigzag direction, respectively. When applying tensile strain along the armchair direction, the $A_{\rm g}^{1}$ mode exhibits a slight redshift, and $B_{\rm 2g}$ and $A_{\rm g}^{2}$ demonstrate a significant blueshift. When applying uniaxial tensile strain along the zigzag direction, all three modes exhibit a redshift.
Subsequently, density functional theory (DFT) calculations have been conducted. In Fig. 1(f), the calculations demonstrate anisotropic band dispersion, consistent with the experimental results. According to the explanation provided by the classical tight-binding model,[17,32] the band dispersion along the Z–L direction is determined by the electron–electron interactions between adjacent phosphorus atom chains along the zigzag direction. In this direction, the distance between atom chains is shorter, and causes stronger electron–electron interactions, leading to pronounced electron hybridization and the formation of band crossings. Meanwhile, the band dispersion along the Z–A direction is determined by the electron–electron interactions between adjacent phosphorus atom chains along the armchair direction. In this direction, the larger distance results in weaker electron hybridization, limiting the formation of band crossings. Furthermore, there are distinctions in the dispersion-$k$ relationships on either side of the high-symmetry point $Z$ in the Brillouin zone. As shown in Fig. S1 of the Supplementary Information, the dispersion-$k_{x}$ relationship along the Z–L (armchair) direction exhibits a quasi-linear relationship [red dashed line in Fig. S1(c)], whereas the dispersion-$k_{y}$ relationship along the Z–A (zigzag) direction follows a quadratic shape [red dashed line in Fig. S1(e)]. These different dispersion-$k$ relationships have been reported in previous research[32,33] and attributed to the variations in the in-plane spin–orbit coupling (SOC) along the Z–L and Z–A directions. Before embarking on the strain-ARPES experiments, the presence of tensile strain in bulk BP was validated by the Raman shift. As depicted in Fig. 2(a), three phonon vibration modes with a relatively strong response exist in the BP crystal: the $A_{\rm g}^{1}$ mode along the out-of-plane direction, the $B_{\rm 2g}$ mode along the zigzag direction, and the $A_{\rm g}^{2}$ mode along the armchair direction. In Figs. 2(b) and 2(c), when applying tensile strain along the armchair direction, the $A_{\rm g}^{1}$ mode exhibits a slight redshift, while the other two modes, $B_{\rm 2g}$ and $A_{\rm g}^{2}$, demonstrate a significant blueshift. Similarly, in Figs. 2(d) and 2(e), when applying tensile strain along the zigzag direction, all three modes exhibit a clear redshift. These results are consistent with previous studies on strain-induced phonon vibrations of thin-layer BP.[27,34-36] Thus, the tensile strain can be applied to the bulk BP along the two special directions. Then, the band structure was investigated by ARPES under the influence of uniaxial stretching, as shown in Fig. S6 of the Supplementary Information.
cpl-41-3-037102-fig3.png
Fig. 3. The VBM shifts of bulk BP regulated by uniaxial tensile strain along the armchair direction. (a)–(f) ARPES mappings, (g) EDC curves, and (h) the VBM shifts. As the tensile strain increases from 0% to 6.25%, the VBM shifts monotonically towards deeper levels, away from the Fermi level.
First, as illustrated in Figs. 3(a)–3(f), a quasi-linear dispersion-$k_{x}$ relationship with a small rounded top remains constant with the tensile strain increasing along the armchair direction. As the strain increases from 0% to 6.25%, the VBM shifts monotonically towards deeper levels, away from the Fermi level (red dashed line with arrow). The corresponding energy distribution curves [EDC, Fig. 3(g)] provide a more intuitive display of the VBM shift, where it shifts from $-210$ meV to $-310$ meV. In Fig. 3(h), the VBM exhibits a downward trend with the increase in tensile strain, and the slope of $-16.2$ meV/% strain is obtained by linear fitting. Similarly, as depicted in Figs. 4(a)–4(i), the quadratic dispersion-$k_{y}$ relationships exist near the top of the valance band (Z–A direction) upon applying the tensile strains. The EDC curves shown in Fig. 4(j) demonstrate that as the strain increases from 0% to 8.33%, the VBM exhibits a significantly monotonic increase from $-210$ meV to $-100$ meV. In Fig. 4(k), a reasonably linear relationship is obtained with a slope of 13.1 meV/% strain.
cpl-41-3-037102-fig4.png
Fig. 4. The VBM shifts of bulk BP regulated by uniaxial tensile strain along the zigzag direction. (a)–(i) ARPES mappings, (j) EDC curves, and (k) the VBM shifts. As the tensile strain increases from 0% to 8.33%, the VBM exhibits a significantly monotonic shift, but this trend gradually moves towards shallower levels.
To further elucidate the observed anisotropic VBM shifts in bulk BP, DFT calculations subjected to uniaxial tensile strain are conducted. The results are presented in Fig. 5(a), which are consistent with the experimental results in Fig. 3(h). To understand the VBM shift, the lattice constants under uniaxial tensile strain are calculated. As shown in Fig. 5(b), the uniaxial tensile strains along the armchair direction (lattice constant $c$) lead to a decrease in the lattice constant $a$ and an increase in the lattice constant $b$. The trend of lattice constant $a$ arises because of the positive Poisson effect in the mutually perpendicular directions within the BP layers. Conversely, the lattice constant $b$ exhibits an unusual expansion because of the negative Poisson effect.[37,38] This enlargement of the interlayer spacing in bulk BP leads to weakening of the interlayer coupling. Based on previous reports,[39-41] the interlayer coupling lowers the bulk bandgap. Thus, the weaker the interlayer coupling, the larger the bandgap, and this shifts the VBM to deeper levels.[32,42] Meanwhile, the effective mass of the electrons and the band velocity near the VBM are derived from the curvature of the bands. The strains along the armchair direction increase the effective mass and reduce the band velocity,[43-45] as shown in Fig. S5 of the Supplementary Information.
cpl-41-3-037102-fig5.png
Fig. 5. The calculated results of the VBM and lattice constant shifts of bulk BP under uniaxial tensile strain along [(a), (b)] the armchair direction and [(c), (d)] the zigzag direction, respectively. Along the armchair direction, the VBM shifts exhibit a consistent monotonic decreasing trend with strain; the increase of the lattice constant $b$ indicates the expansion of the interlayer spacing in bulk BP. Along the zigzag direction, the VBM shift trend is completely opposite to that along the armchair direction; the decrease of the lattice constant $b$ indicates the compression of the interlayer spacing in bulk BP.
Similarly, the calculated results along the zigzag direction, as shown in Fig. 5(c), exhibit VBM behavior opposite to that along the armchair direction, which is also consistent with Fig. 4(k). The tensile-strain-induced lattice constant shifts are presented in Fig. 5(d). When the tensile strains are applied along the zigzag direction (lattice constant $a$), the lattice constants $b$ and $c$ decrease because of the positive Poisson effect. The decrease in the interlayer spacing enhances the interlayer coupling and lowers the bandgap, and the VBM shifts towards shallower levels. As depicted in Fig. S5, with the increase in strain along the zigzag direction, it results in a smaller effective mass and a faster band velocity. This work clearly and comprehensively reveals the band structure of bulk BP along the Z–L (armchair) and Z–A (zigzag) directions. In the Z–L direction, strong intralayer electron–electron interactions result in a bowknot-like band structure with a quasi-linear dispersion-$k_{x}$ relationship. Conversely, the dispersion-$k_{y}$ relationships exhibit a quadratic shape without band crossings in the Z–A direction due to weaker intralayer electron–electron interactions. Subsequently, anisotropic band evolutions are investigated under uniaxial tensile strains along the armchair and zigzag directions. The results indicate that increasing tensile strain along the armchair direction (from 0% to 6.25%) expands interlayer spacing, weakens interlayer coupling, and shifts the VBM to deeper levels, with a slope of $-16.2$ meV/% strain. When the tensile strain along the zigzag direction increases from 0% to 8.33%, the interlayer spacing decreases, enhances the interlayer coupling, and shifts the VBM to shallower levels, with a slope of 13.1 meV/% strain. These results elucidate the relationship between the interlayer spacing and the VBM in bulk BP, offering valuable insights for band engineering and potential applications. Acknowledgments. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12104216, 12241403, and 61974061), the National Key R&D Program of China (Grant No. 2021YFB3601600), and the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20140054). References High-mobility transport anisotropy and linear dichroism in few-layer black phosphorusUnconventional conductivity increase in multilayer black phosphorusInfrared fingerprints of few-layer black phosphorusBlack phosphorus field-effect transistorsStrain-Modulated Bandgap and Piezo-Resistive Effect in Black Phosphorus Field-Effect TransistorsPerformance Improvement in Hydrogenated Few-Layer Black Phosphorus Field-Effect TransistorsDetermination of layer-dependent exciton binding energies in few-layer black phosphorusBand Gap Renormalization, Carrier Multiplication, and Stark Broadening in Photoexcited Black PhosphorusStrain-engineered widely tunable perfect absorption angle in black phosphorus from first principlesRecent progress in black phosphorus and black-phosphorus-analogue materials: properties, synthesis and applicationsSemiconducting black phosphorus: synthesis, transport properties and electronic applicationsTunable optical properties of multilayer black phosphorus thin filmsElectrically Tunable Energy Bandgap in Dual-Gated Ultra-Thin Black Phosphorus Field Effect TransistorsSuperior mechanical flexibility of phosphorene and few-layer black phosphorusElastic properties of suspended black phosphorus nanosheetsScaling laws of band gaps of phosphorene nanoribbons: A tight-binding calculationElectronic Structure of Black Phosphorus in Tight Binding ApproachElectronic Structure of Black Phosphorus in Self-Consistent Pseudopotential ApproachTunable Electrical Performance of Few-Layered Black Phosphorus by StrainExciton in phosphorene: Strain, impurity, thickness, and heterostructureLayer-controlled band gap and anisotropic excitons in few-layer black phosphorusLayer-Dependent Pressure Effect on the Electronic Structure of 2D Black PhosphorusStrain Effects on Properties of Phosphorene and Phosphorene Nanoribbons: a DFT and Tight Binding StudyPressure-Induced Electronic Transition in Black PhosphorusStrain-Engineering the Anisotropic Electrical Conductance of Few-Layer Black PhosphorusAnisotropic Electron–Phonon Interactions in Angle-Resolved Raman Study of Strained Black PhosphorusGiant Anisotropic Raman Response of Encapsulated Ultrathin Black Phosphorus by Uniaxial StrainResolving the In‐Plane Anisotropic Properties of Black PhosphorusAnisotropic Electron-Photon and Electron-Phonon Interactions in Black PhosphorusProlonged photo-carriers generated in a massive-and-anisotropic Dirac materialNon-equilibrium band broadening, gap renormalization and band inversion in black phosphorusObservation of tunable band gap and anisotropic Dirac semimetal state in black phosphorusElectronic structure of bulk and multilayer black phosphorus under strain: a minimal model studyIn‐Plane Uniaxial Strain in Black Phosphorus Enables the Identification of Crystalline OrientationLattice Vibration and Raman Scattering in Anisotropic Black Phosphorus CrystalsLattice vibrational modes and Raman scattering spectra of strained phosphoreneNegative poisson’s ratio in single-layer black phosphorusAuxetic Black Phosphorus: A 2D Material with Negative Poisson’s RatioDirect observation of the layer-dependent electronic structure in phosphoreneEffect of van der Waals interactions on the structural and elastic properties of black phosphorusStrain-tunable van der Waals interactions in few-layer black phosphorusStrain-controlled fundamental gap and structure of bulk black phosphorusPerformance analysis of uniaxially strained monolayer black phosphorus and blue phosphorus n-MOSFET and p-MOSFETAnisotropic basic electronic properties of strained black phosphoreneElectronic structure and optic absorption of phosphorene under strain
[1] Qiao J S, Kong X H, Hu Z X, Yang F, and Ji W 2014 Nat. Commun. 5 4475
[2] Koleśnik-Gray M, Meingast L, Siebert M, Unbehaun T, Huf T, Ellrott G, Abellán G, Wild S, Lloret V, Mundloch U, Schwarz J, Niebauer M, Szabo M, Rommel M, Hutzler A, Hauke F, Hirsch A, and Krstić V 2023 npj 2D Mater. Appl. 7 21
[3] Zhang G W, Huang S Y, Chaves A, Song C Y, ÖzÇelik V O, Low T, and Yan H G 2017 Nat. Commun. 8 14071
[4] Li L K, Yu Y J, Ye G J, Ge Q Q, Ou X D, Wu H, Feng D L, Chen X H, and Zhang Y B 2014 Nat. Nanotechnol. 9 372
[5] Zhang Z C, Li L K, Horng J, Wang N Z, Yang F Y, Yu Y J, Zhang Y, Chen G R, Watanabe K, Taniguchi T, Chen X H, Wang F, and Zhang Y B 2017 Nano Lett. 17 6097
[6] Zheng H M, Sun S M, Liu H, Huan Y W, Yang J G, Zhu B, Liu W J, and Ding S J 2018 Chin. Phys. Lett. 35 127302
[7] Zhang G W, Chaves A, Huang S Y, Wang F J, Xing Q X, Low T, and Yan H G 2018 Sci. Adv. 4 eaap9977
[8] Chen Z S, Dong J W, Papalazarou E, Marsi M, Giorgetti C, Zhang Z L, Tian B, Rueff J P, Taleb-Ibrahimi A, and Perfetti L 2019 Nano Lett. 19 488
[9] Alidoust M, Halterman K, Pan D, Willatzen M, and Akola J 2020 Phys. Rev. B 102 115307
[10] Xu Y J, Shi Z, Shi X Y, Zhang K, and Zhang H 2019 Nanoscale 11 14491
[11] Liu H, Du Y C, Deng Y X, and Ye P D 2015 Chem. Soc. Rev. 44 2732
[12] Low T, Rodin A S, Carvalho A, Jiang Y J, Wang H, Xia F, and Castro A H 2014 Phys. Rev. B 90 075434
[13] Yan S L, Xie Z J, Chen J H, Taniguchi T, and Watanabe K 2017 Chin. Phys. Lett. 34 047304
[14] Wei Q and Peng X H 2014 Appl. Phys. Lett. 104 251915
[15] Wang J Y, Li Y, Zhan Z Y, Li T, Zhen L, and Xu C Y 2016 Appl. Phys. Lett. 108 013104
[16] Taghizadeh S E, Zare M H, and Fazileh F 2015 Phys. Rev. B 91 085409
[17] Takao Y and Morita A 1981 J. Phys. Soc. Jpn. 50 3362
[18] Asahina H, Shindo K, and Morita A 1982 J. Phys. Soc. Jpn. 51 1193
[19] Wu H R, Liu X F, Yin J, Zhou J X, and Guo W L 2016 Small 38 5276
[20] Arra S, Babar R, and Kabir M 2019 Phys. Rev. B 99 045432
[21] Tran V, Soklaski R, Liang Y, and Yang L 2014 Phys. Rev. B 89 235319
[22] Huang S Y, Lu Y, Wang F J, Lei Y C, Song C Y, Zhang J S, Xing Q X, Wang C, Xie Y G, Mu L, Zhang G W, Yan H, Chen B, and Yan H G 2021 Phys. Rev. Lett. 127 186401
[23] Zhang R Y, Zheng J M, and Jiang Z Y 2018 Chin. Phys. Lett. 35 017302
[24] Xiang Z J, Ye G J, Shang C, Lei B, Wang N Z, Yang K S, Liu D Y, Meng F B, Luo X G, Zou L J, Sun Z, Zhang Y, and Chen X H 2015 Phys. Rev. Lett. 115 186403
[25] Fei R X and Yang L 2014 Nano Lett. 14 2884
[26] Zhu W N, Liang L B, Roberts R H, Lin J F, and Akinwande D 2018 ACS Nano 12 12512
[27] Li Y Y, Hu Z X, Lin S H, Lai S K, Ji W, and Lau S P 2017 Adv. Funct. Mater. 27 1600986
[28] Liu X L, Ryder C R, Wells S A, and Hersam M C 2017 Small Methods 1 1700143
[29] Ling X, Huang S X, Hasdeo E H, Liang L B, Parkin W M, Tatsumi Y, Nugraha A R T, Puretzky A A, Das P M, Sumpter B G, Geohegan D B, Kongv J, Saito R, Drndic M, Meunier V, and Dresselhaus M S 2016 Nano Lett. 16 2260
[30] Nurmamat M, Ishida Y, Yori R, Sumida K, Zhu S, Nakatake M, Ueda Y, Taniguchi M, Shin S, Akahama Y, and Kimura A 2018 Sci. Rep. 8 9073
[31] Hedayat H, Ceraso A, Soavi G, Akhavan S, Cadore A, Dallera C, Cerullo G, Ferrari A C, and Carpene E 2021 2D Mater. 8 025020
[32] Kim J, Baik S S, Ryu S H, Sohn Y, Park S, Park Y G, Denlinger J, Yi Y, Choi H J, and Kim K S 2015 Science 349 723
[33] Yang M and Liu H 2020 Phys. Scr. 95 035805
[34] Zhang Q, Mao N N, Wu J X, Tong L M, Zhang J, and Liu Z R 2017 Small 13 1700466
[35] Mao N N, Zhang S Q, Wu J X, Zhang J, and Tong L M 2018 Small Methods 2 1700409
[36] Fei R and Yang L 2014 Appl. Phys. Lett. 105 083120
[37] Jiang J W and Park H S 2014 Nat. Commun. 5 4727
[38] Du Y C, Maassen J, Wu W R, Luo Z, Xu X, and Ye P D 2016 Nano Lett. 16 6701
[39] Li L K, Kim J, Jin C H, Ye G J, Qiu D Y, Jornada F H, Shi Z W, Chen L, Zhang Z C, Yang F Y, Watanabe K, Taniguchi T, Ren W C, Louie S G, Chen X H, Zhang Z, and Wang F 2017 Nat. Nanotechnol. 12 21
[40] Appalakondaiah S, Vaitheeswaran G, Lebègue S, Christensen N E, and Svane A 2012 Phys. Rev. B 86 035105
[41] Huang S Y, Zhang G W, Fan F R, Song C Y, Wang F J, Xing Q X, Wang C, Wu H, and Yan H G 2019 Nat. Commun. 10 2447
[42] Guan J, Song W, Yang L, and Tomanek D 2016 Phys. Rev. B 94 045414
[43] Banerjee L, Mukhopadhyay A, Sengupta A, and Rahaman H 2016 J. Comput. Electron. 15 919
[44] Yarmohammadi M, Mortezaei M, and Mirabbaszadeh K 2020 Physica E 124 114323
[45] Duan H J, Yang M, and Wang R Q 2016 Physica E 81 177