Anisotropic exchange in GdGa

Anisotropic exchange in GdGa

Journal of Alloys and Compounds 827 (2020) 154119 Contents lists available at ScienceDirect Journal of Alloys and Compounds journal homepage: http:/...

558KB Sizes 3 Downloads 47 Views

Journal of Alloys and Compounds 827 (2020) 154119

Contents lists available at ScienceDirect

Journal of Alloys and Compounds journal homepage: http://www.elsevier.com/locate/jalcom

Anisotropic exchange in GdGa  brega , P.O. Ribeiro , B.P. Alho , P.J. von Ranke V.S.R. de Sousa *, E.P. No Instituto de Física Armando Dias Tavares, Universidade do Estado do Rio de Janeiro (UERJ), 20550-013, Rio de Janeiro - RJ, Brazil

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 October 2019 Received in revised form 28 January 2020 Accepted 29 January 2020 Available online 1 February 2020

In this work we discuss the non-collinear ground-state reported for GdGa on the basis of a model Hamiltonian considering anisotropic exchange interactions. We show that a competition between intra and inter-sublattice exchange interactions can lead to a canted structure, which competes with a collinear (ferro or antiferromagnetic) order. The mean-field thermodynamic analysis of the model for S ¼ 7=2 shows a good agreement between calculated and reported experimental data of the canting angle and isothermal entropy change for GdGa. © 2020 Elsevier B.V. All rights reserved.

Keywords: Spin Hamiltonians Anisotropic exchange Magnetocaloric effect

1. Introduction GdGa crystallizes in the orthorhombic CrB-type structure, space group Cmcm (#63), in which Gd and Ga occupy 4c site [1e4]. It orders ferromagnetically around 200 K and undergoes a spinreorientation transition at a lower temperature. The ferromagnetic order is collinear between Tc and TSR, and for temperatures below TSR it becomes non-collinear, with Gd 4c site being split into two non-equivalent ones. This spin reorientation has been studied €ssbauer spectroscopy and neutron experimentally with both Mo powder diffraction [4,5]. Since Gd3þ have a spherical symmetrical charge distribution, resulting from the half-filled 4f-shell configuration, we can rule out the influence of the crystal field interaction on the spin reorientation in GdGa. One may speculate about the influence of the classical dipole-dipole interaction [6,7], which leads to anisotropic effects, but in the present case the reported reorientation temperatures are too high to be explained only by this interaction. The antisymmetric exchange [8,9], which tends to orient two interacting spins perpendicularly, destabilizing collinear ferromagnetic or antiferromagnetic structures, and can cause a canting of the spins in different sites [10], may be ruled out on a symmetry basis: since the middle of the bonds between Gd-atoms in GdGa are inversion centers the DzyaloshinskiieMoriya interaction is null [11]. As pointed out by Delyagin et al. [5], this reorientation might be a consequence of anisotropic exchange interactions, which may be

* Corresponding author. E-mail address: [email protected] (V.S.R. de Sousa). https://doi.org/10.1016/j.jallcom.2020.154119 0925-8388/© 2020 Elsevier B.V. All rights reserved.

related with an anisotropic GdeGa hybridization mediating the indirect GdeGd interaction. Here we report theoretical results, that account for the major experimental features of GdGa, based on a microscopic model Hamiltonian that considers anisotropic exchange interactions between Gd ions. The model Hamiltonian is a two-sublattice anisotropic Heisenberg model in the presence of a magnetic field. In order to apply the model to GdGa we consider its magnetic structure as composed of alternating ferromagnetic planes as speculated by Leithe-Jasper and Hiebl [12]. The exchange parameters were obtained semi-empirically following a detailed analysis of the system behavior near the phase transitions. The calculated spin canting angle and isothermal entropy change show a good agreement with the reported experimental data for GdGa.

2. Formalism We consider the following anisotropic Heisenberg model in the presence of a magnetic field:

H ¼

 ! X X ⊥ z z k  x x ! y y jij Si Sj þ jij Si Sj þ Si Sj  B 0 , gi mB S i : ij

(1)

i

! a z is the where jij (a ¼ ⊥; k) are the exchange parameters, B 0 ¼ B0 b  factor of the i-th spin and mB the applied magnetic field, gi the Lande Bohr magneton. We also assume a magnetic structure composed of alternating ferromagnetic layers in which the spins are either parallel to a global z-axis or canted from this axis by an angle q. This arrangement is sketched in Fig. 1.

2

V.S.R. de Sousa et al. / Journal of Alloys and Compounds 827 (2020) 154119

(i) j2 > j1 favors FM and j2 < j1 favors a canted (non-collinear) magnetic structure. Conversely, when j2 is negative the FM phase is always unstable and for the case j1 > 0 one may observe a competition between AFM (favored if jj2 j > j1 ) and NCM (favored if jj2 j < j1 ) states. The classical ground state phase diagram is depicted in Fig. 2. The transition between the different magnetic states is marked by the reversible lines j2 ¼ j1 (FM % NCM), j2 ¼ j1 (AFM % NCM) and j2 ¼ 0 (FM % AFM, with j1 < 0Þ). There is a critical magnetic field (Bc ), applied along z, for which the FM state becomes more favorable than the NCM one, even if j1 > j2 . One can show from (3) that this field will be given by

Bc ¼

2ðj2  j1 Þ S: g mB

(4)

This classical picture illustrates that an anisotropic exchange interaction may be responsible for the non-collinear magnetic state observed in GdGa. Fig. 1. Sketch of the layered magnetic structure considered for the model Hamiltonian.

For the sake of simplicity, the exchange interaction parameters k between spins in both sublattices are taken as: j11 ¼ j1 , j⊥ 11 ¼ 0; k k ⊥ j22 ¼ 0, j⊥ ¼ j’ ; j ¼ j , j ¼ 0. This implies that the intra-layer 1 12 2 12 22 exchange between spins in the layer composed of spins parallel to z is j1 ; the intra-layer exchange between canted spins is along an axis perpendicular to z; and the inter-layer exchange j2 is along z. For convenience, we name the spins in the layer composed of spins ! ! parallel to z as S 1 , and the canted ones as S 2 , and we also limit the calculations to the X Z plane. Considering the above set of exchange parameters, we may decompose relation (1) in two-sublattice Hamiltonians, i.e., H ¼ H 1 þ H 2 , so that

H

1

¼  j1

X

Sz1i Sz1j  j2

X

CijD

H

2

¼  j’1

X

Sz1i Sz2j  mB B0

X

CijD

Sx2i Sx2j  j2

CijD

X CijD

g1i Sz1i ;

(2a)

g2i Sz2i :

(2b)

i

Sz2i Sz1j  mB B0

X

2.2. Mean field analysis for S ¼ 1=2 Here we discuss the above model for an ensemble of spins with S ¼ 1=2 in the mean field approximation. This allows us to obtain equations of state that, in turn, can be evaluated at some critical conditions in order to establish possible phase transitions. In the mean-field approximation, equation (2) are rewritten in the form of an effective field acting on each sublattice, so that, H 1 ¼ b1z Sz1 and H 2 ¼  b2z Sz2  b2x Sx2 , where b1z ¼ b0 þ 2z1 j1 CSz1 D þ 2z2 j2 CSz2 D, b2z ¼ b0 þ 2z2 j2 CSz1 D, b2x ¼ 2z1 j1 CSx2 D, b0 ¼ g mB B0 . The eigenvalues of H ð1Þ

Sb1z , ε2

Once we have obtained the eigenvalues, we can proceed by calculating the mean thermodynamic values of the operators Sz1 , Sx2 and Sz2 using the relation:

CSkl D ¼

i

In the above relations the summation is restricted to nearestneighbors (NN). Note that in the present form, relations (2) could be used to study a layered structure of unequal moments. Here we make the further assumption that S1 ¼ S2 ¼ S (g1 ¼g2 ¼ g) and j’1 ¼ j1 , this last one is justifiable since we are considering a magnetic structure of equal magnetic moments.

ð2Þ

¼ þ Sb1z , ε1

ð1Þ

and H 2 will be given by ε1 ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ ¼  S b22x þ b22z , ε2 ¼ þ S b22x þ b22z . 1

1X Zl l



! ðlÞ ðlÞ vεi ebεi ; vblk

(5)

P where l ¼ 1,2, k ¼ x,z,. Zl ¼ ebεi ðlÞ is the canonical partition i function, and b ¼ 1=kB T (kB is the Boltzmann constant and T the temperature). After some algebra one gets the following results

2.1. Classical ground state In a classical picture, the spin operators in both sublattices will ! ! be given by: S 1 ¼ Sz1 b z ¼ Sb z ; S 2 ¼ Sx2 b z , Sx2 ¼ Ssinq and x þ Sz2 b z S2 ¼ Scosq. The total energy E ¼ E1 þ E2 is therefore

  E ¼  j1 S2 1 þ sin2 q  2j2 S2 cosq  g mB B0 Sð1 þ cosqÞ:

(3)

By minimizing this energy with respect to q one gets the following critical angles in the absence of an external magnetic field: q ¼ 0, q ¼ p and q ¼ cos1 ðj2 =j1 Þ. Note that this gives two collinear solutions, one corresponding to ferromagnetism (q ¼ 0) and the other to antiferromagnetism (q ¼ p); and a non-collinear solution that depends on the ratio j2 =j1 . For the case in which both j2 and j1 are positive the antiferromagnetic (AFM) solution is always unstable, whereas there will be a competition between ferromagnetic (FM) and non-collinear magnetic (NCM) states since:

Fig. 2. Classical ground state magnetic phase diagram. The competition between the different magnetic structures is depicted as a function of j1 and j2 .

V.S.R. de Sousa et al. / Journal of Alloys and Compounds 827 (2020) 154119

CSz1 D ¼ S tan hðbSb1z Þ;

(6a)

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  b2x ffi tan h bS b22x þ b22z ; CSx2 D ¼ S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b22x þ b22z

(6b)

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  b2z ffi tan h bS b22x þ b22z : CSz2 D ¼ S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b22x þ b22z

(6c)

(7)

where a ¼ s1 z2 j2 =z1 j1 and s1 ¼ CS1z D is the expectation value of S1Z evaluated from Eq. (6a) at TSR. This implies that equation (7) is calculated self-consistently, in a similar fashion as the critical temperature of the Oguchi method [13,14]. One may obtain the el) temperature from relations ð6aÞ and (6c) in the limit Curie (or Ne of small arguments, so that

4S2 a z2 j2 TCðNÞ ¼  : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  kB  1 þ 1 þ 4a2

In order to obtain the spin reorientation temperature for any Svalue, we consider the Hamiltonian

H ¼  b2z Sz2  b2x Sx2 ;

(9)

(10)

as composed of an unperturbed

H

0

¼  b2z Sz2

(11)

and a perturbed

H ’ ¼  b2x Sx2 ¼ 

b2x þ ðS þ S Þ 2

(12)

contribution. This is well justified since close to TSR the perpendicular magnetization (and, henceforth, the effective field along x) ð0Þ approaches zero. The unperturbed eigenenergies εm ¼ b2z m  are ð0Þ easily obtained, with the corresponding eigenvectors εm D ¼ S; mD (m ¼  S; …; S). The first order correction to the eigenenergies due ð1Þ ð0Þ ð0Þ to the perturbation εm ¼ Cεm jH 1 jεm D is null, and to second order the correction will be given by ð2Þ

εm ¼ 

b22x m: 2b2z

(13)

The mean thermodynamic value of S2x can then be obtained from relation

(8)

One may also consider the limit T/0. In this case, CS1z D ¼ S, and we can define the ratio between CSx2 D and CSz2 D as the tangent of the canted angle. After some algebra, one gets

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   z1 =z2 2 tanq ¼  1: j2 =j1

In the above relation the canting angle depends on the coordination number as shown in Fig. 3.

2.3. Critical temperatures for S > 1=2

Relations (6) form a group of magnetic equations of state that can be used to evaluate the dependence of the magnetization as a function of temperature and magnetic field. Here we are mainly interested in the behavior around the spin reorientation temperature (TSR ) and the Curie temperature (TC ). TSR marks the transition from the noncollinear to a collinear magnetic state (either FM if j2 > 0, or AFM if j2 < 0). At this temperature the magnetization of each sublattice becomes parallel to the z-axis, which implies that CSx2 D /0 at TSR . Using this condition in (6b) one finds that (for B0 ¼ 0) the reorientation temperature has the following dependence on the exchange parameters

4 s S2 z j TSR ¼  1  2 2 ; kB ln 1þa 1a

3

CS2x D ¼

S  X

1 Z ð0Þ

m¼S



ð2Þ  ð0Þ vεm ebεm ; vb2x

(14)

P ð0Þ where Z ð0Þ ¼ Sm¼S ebεm . One can show that the above equation can be cast in the form

CS2x D ¼ S

b2x B ðbSb2z Þ; b2z S

(15)

where BS ðxÞ is the Brillouin function. This result is valid close to the spin reorientation temperature. We can therefore substitute b2x and b2z above and consider T ¼ TSR to get the following equation

TSR ¼

2s1 S2 z2 j2 : kB B1 S ðaÞ

(16)

We may now consider the following expansion for the inverse Brillouin function [15].

B1 S ðaÞ ¼

kðaÞ 2

ln

  1þa ; 1a

(17)

where

kðaÞ ¼ Fig. 3. Canting angle as a function of the ratio j2 =j1 for several coordinations z1 = z2 ¼ 2 (simple cubic and simple tetragonal lattices), 1/2 (face-centered and body-centered cubic lattices), 3 (hexagonal lattice). q ¼ 90+ implies non-interacting sublattices and q ¼ 0 ferromagnetism.

15  11ð1  εÞð1 þ 2εÞa2 5 þ 10ε  ð1  εÞ½5 þ 11εð1 þ 2εÞa2

(18)

with ε ¼ 1=2S. Therefore the spin reorientation temperature for any S value will be given by

4

V.S.R. de Sousa et al. / Journal of Alloys and Compounds 827 (2020) 154119

TSR ¼

4s1 S2 z j   2 2; kB 1þa kðaÞln 1a

(19)

which differs from (7) by a factor kðaÞ in the denominator. Note that, if ε ¼ 1/kðaÞ ¼ 1, and we recover the result obtained for S ¼ 1= 2. One may also show that the equation for TCðNÞ is generalized by replacing S2 by SðS þ1Þ=3 in (8), therefore

TCðNÞ ¼

4SðS þ 1Þ a  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi z2 j2 : 3kB  1 þ 1 þ 4a2

(20)

Relations (7) and (8), or (19) and (20), allows one to obtain for some experimentally determined spin reorientation and critical temperatures the values of j2 and j1 . These parameters can be used as input values in the model and one may evaluate the equations of state for the magnetization, and also obtain other thermodynamic quantities (e.g., specific heat) as a function of temperature and magnetic field. With these semi-empirical exchange parameters one can numerically estimate the canting angle using

q ¼ tan1

  CS2x D CS2z D

(21)

We point out that the critical magnetic field, at which the NCM state becomes less stable than the FM one, will also be given by (4) with a slight modification j1 /z1 j1 and j2 /z2 j2 . 3. Results and discussions Here we consider the model described in section 2 to study the spin reorientation transition and the thermomagnetic behavior of the compound GdGa, which shows a non-collinear to ferromag factor and spin angular netic phase transition. We take the Lande momentum given by Hund’s rules: g ¼ 2 and S ¼ 7= 2. From the reported crystallographic structure of GdGa, gadolinium atoms have two first neighbors and 4 s neighbors, therefore z1 ¼ 2 and z2 ¼ 4. We have calculated the isothermal entropy change by the standard procedure described, for instance, in Ref. [16]. The magnetic entropy of each sublattice is obtained from the relation Smag ¼ NkB ðlnZl þ bEl Þ, where El ¼ CH l D represents the mean thermodynamic value of the l-th Hamiltonian in (2). In order to compare the calculations of the isothermal entropy change with experimental data we adopt a procedure in which the applied magnetic field direction is varied over a unit sphere, in an averaging procedure used to simulate a polycrystalline sample [17]. Table 1 lists some reported experimental critical temperatures for GdGa as well as the experimental relative canting angle between Gd atoms in the non-collinear phase (qexp ). We also list the corresponding exchange parameters obtained from equations (19) and (20) and the canting angle calculated from (21). One can note a good agreement between calculated and experimental canting angles, showing that in fact the mechanism responsible for the

Table 1 Experimental Curie temperature (TC), spin reorientation temperature (TSR) and relative angle between Gd moments (qexp ), and calculated intra-sublattice exchange parameter (j1 ), inter-sublattices exchange parameter (j2 ) and canting angle (qcalc ). TC

TSR

qexp ( )

Ref.

j1 (meV)

j2 (meV)

qcalc ( )

190 183

68 100

38 45

[4] [5]

0.514 0.533

0.226 0.199

28 42

non-collinear magnetic state in GdGa is an anisotropic exchange. Fig. 4 shows the calculated zero-field reduced sublattices magnetization (sk ¼CSk D=S) as a function of temperature using the exchange parameters j1 ¼ 0:533 meV and j2 ¼ 0:199 meV. At very low temperatures, s1z is fully saturated and the components of the canted sublattice have the values s2z ¼ cosqcalc and s2x ¼ sinqcalc . As temperature increases the magnetization of each sublattice decreases, with the transverse magnetization (parallel to x) going to zero at TSR. Above the reorientation temperature, the system evolves as two coupled ferromagnetic sublattices with different magnetic moments. The behavior of the curves close to the critical temperatures is unaffected by short-range correlations since we are considering a mean-field approximation, therefore the order parameter goes to zero above these temperatures. The inset in Fig. 4 shows the variation of the canting angle with temperature. At TSR both sublattices become parallel to the z axis and, therefore, q ¼ 0 above this temperature. The smooth change of the canting angle as a function of temperature shows that the spin reorientation transition is of second-order, which is in accordance with an Arrott plot analysis [18,19] performed by Zhang et al. [3], which also studied the magnetocaloric properties of GdGa, reporting a maximum value of 4.81 J/kg.K for the isothermal entropy change under a magnetic field change from 0 to 5 T. The calculated entropy change shows a good agreement with the experimental data as can be seen in Fig. 5. The secondary peak presented in the theoretical curves are due to the spin reorientation transition, which has a similar impact on the isothermal entropy change as observed in the compounds Ho2In Ref. [20], ErGa [21], TbZn [22], HoZn [23,24], i.e., a large relative cooling power (RCP) over a wide temperature range [3]. The main difference being that in these compounds the spin reorientation appears as a competition between crystal field and exchange interactions, meanwhile in GdGa it is due only to the anisotropic exchange mechanism. 4. Summary and outlook We have shown that the spin reorientation transition and the large refrigeration capacity reported for GdGa may be explained by an anisotropic two-coupled sublattices Heisenberg-like model Hamiltonian. Our mean-field results are in good agreement with

Fig. 4. Zero field reduced magnetization sublattices as a function of temperature calculated using j1 ¼ 0:533 meV and j2 ¼ 0:199 meV. The vertical dashed lines mark the reported experimental TSR and the shaded area the range of reported experimental Curie temperatures for GdGa. The inset shows the evolution of the canting angle as a function of temperature.

V.S.R. de Sousa et al. / Journal of Alloys and Compounds 827 (2020) 154119

5

CRediT authorship contribution statement V.S.R. de Sousa: Conceptualization, Methodology, Software, Formal analysis, Writing - original draft, Funding acquisition,  brega: Validation, Writing - review Project administration. E.P. No & editing, Funding acquisition. P.O. Ribeiro: Validation, Writing review & editing, Funding acquisition. B.P. Alho: Validation, Writing - review & editing, Funding acquisition. P.J. von Ranke: Methodology, Validation, Writing - review & editing, Funding acquisition, Project administration. Acknowledgements The authors thank Brazilian agencies FAPERJ and CNPq for financial support. This study was financed in part by the Coordenaç~ ao de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. Fig. 5. GdGa isothermal entropy change as a function of temperature for several magnetic field changes. Full circles and solid lines represent, respectively, experimental data and calculations.

the experimental data for both canting angles and isothermal entropy change. The canting angle is dependent on the ratio between intra and intersublattice exchanges, and the anisotropic exchange driven spin reorientation transition is responsible for the large observed RCP reported for GdGa [3]. There is a previous theoretical study that reports density functional theory calculated exchange interactions [25] for GdGa [26], which cannot be directly compared to the exchange parameters adopted here because the authors have considered only a simple collinear ferromagnetic structure in their study, not considering the influence of the spin-orbit interaction in their calculations. Therefore, it would be interesting if further investigations were performed in order to describe the non-collinear magnetic ground state of GdGa from first-principles, in order to clarify the possible anisotropic hybridization between Gd and Ga atoms in this compound. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1] N.C. Baenziger, J.L. Moriarty Jr., Acta Crystallogr. 14 (1961) 946. [2] K.H.J. Buschow, W.W.v. d. Hoogenhof, J. Less. Comm. Met. 45 (1976) 309. [3] J.Y. Zhang, J. Luo, J.B. Li, J.K. Liang, Y.C. Wang, L.N. Ji, Y.H. Liu, G.H. Rao, J. Alloys Compd. 469 (2009) 15. ~ oz Pe rez, [4] R.A. Susilo, J.M. Cadogan, D.H. Ryan, N.R. Lee-Hone, R. Cobas, S. Mun Hyperfine Interact. 226 (2014) 257. [5] N. Delyagin, V. Krylov, I. Rozantsev, J. Magn. Magn Mater. 308 (2007) 74. [6] J.M. Luttinger, L. Tisza, Phys. Rev. 70 (1946) 954. [7] D.C. Johnston, Phys. Rev. B 93 (2016), 014421. [8] I. Dzyaloshinsky, J. Phys. Chem. Solid. 4 (1958) 241. [9] T. Moriya, Phys. Rev. Lett. 4 (1960) 228. [10] S. Blügel, G. Bihlmayer, Magnetism of low-dimensional systems: Theory, in: €hnle, S. Maekawa, I. Zutic (Eds.), Handbook of H. Kronmüller, S. Parkin, M. Fa Magnetism and Advanced Magnetic Materials, Wiley-Interscience, 2007. pieux, C. Lacroix, J. Magn. Magn Mater. 182 (1998) 341. [11] A. Cre [12] A. Leithe-Jasper, K. Hiebl, Phys. Stat. Sol. 155 (1996) 223. [13] T. Oguchi, Prog. Theor. Phys. 13 (1955) 148. [14] J.S. Smart, Effective Field Theories of Magnetism, Saunders, London, 1966. €ger, J. Non-Newtonian Fluid Mech. 223 (2015) 77. [15] M. Kro [16] N.A. de Oliveira, P.J. von Ranke, Phys. Rep. 489 (2010) 89. [17] B.P. Alho, A.M.G. Carvalho, P.J. von Ranke, J. Magn. Magn Mater. 324 (2012) 210. [18] A. Arrott, Phys. Rev. 108 (1957) 1394. [19] I. Yeung, R.M. Roshko, G. Williams, Phys. Rev. B 34 (1986) 3456. [20] Q. Zhang, J.H. Cho, B. Li, W.J. Hu, Z.D. Zhang, Appl. Phys. Lett. 94 (2009), 182501. [21] J. Chen, B.G. Shen, Q.Y. Dong, F.X. Hu, J.R. Sun, Appl. Phys. Lett. 95 (2009) 132504. [22] V.S.R. de Sousa, E.J.R. Plaza, P.J. von Ranke, J. Appl. Phys. 107 (2010), 103928. [23] V.S.R. de Sousa, E.J.R. Plaza, F.C.G. Gandra, J. Appl. Phys. 109 (2011), 063904. €ttgen, S. Zhou, J. Alloys Compd. 643 (2015) 147. [24] L. Li, Y. Yuan, Y. Zhang, R. Po [25] M.I. Katsnelson, A.I. Lichtenstein, Phys. Rev. B 61 (2000) 8906. [26] X.B. Liu, Z. Altounian, Physica B 406 (2011) 710.