You are an expert at summarizing long articles. Proceed to summarize the following text: 
recently unconventional electronic properties of mono - crystalline graphenes attracts much attention motivated by experimental fabrication,@xcite although they were already the subject of theoretically study prior to the fabrication.@xcite multilayer films which contain more than two layers can also be synthesized , and various phenomena depending on the layer number have been reported.@xcite in this paper we present a theoretical study on the orbital magnetism in multilayer graphenes .    the electronic structure of the monolayer graphene is quite different from conventional metals , because the conduction and valence bands touch at @xmath1 and @xmath2 points in the brillouin zone , around which the dispersion becomes linear like a relativistic particle . in multi - layer graphenes , the interlayer coupling makes a complex structure around the band touching . the electronic properties of graphene bilayer were theoretically studied for the band structure @xcite and the transport properties.@xcite for few - layered graphenes of more than two stacks , the electronic structure is investigated theoretically in a @xmath3 approximation,@xcite a density functional calculation,@xcite and a tight - binding model.@xcite on the experimental side , the band structures of graphenes from one to four layers were recently measured using angle - resolved photoemission spectroscopy.@xcite    the orbital magnetism in graphene - based systems was first studied for a monolayer as a simple model to explain the large diamagnetism of graphite.@xcite it was found that the susceptibility becomes highly diamagnetic at @xmath4 ( band touching point ) even though the density of states vanishes there . the calculation was extended to graphite @xcite and to few - layered graphenes as a model of graphite intercalation compounds.@xcite the fermi surface of the graphite is known to be trigonally warped around the band touching point @xcite and the effect of the warping on magnetization was discussed within the perturbational approach.@xcite recently , the disorder effects on the magnetic oscillation @xcite and on the susceptibility @xcite were studied for the monolayer graphene .    here we present a systematic study on the orbital magnetism for multilayer graphenes with arbitrary layer numbers in the effective mass approximation . we show that the hamiltonian of a multilayer graphene can be decomposed into those equivalent to monolayer or bilayer , which allows us to study the dependence of the susceptibility on layer numbers . we take the trigonal warping effect into the calculation and show that the fine structure around zero energy gives rise to singular magnetic properties . we introduce the model hamiltonian and its decomposition into subsystems in sec . [ sec_form ] , and present the calculation of the magnetization in sec . [ sec_mag ] . the discussion and summary are given in sec . [ sec_disc ] . we consider a multilayer graphene composed of @xmath5 layers of a carbon hexagonal network , which are arranged in the ab ( bernal ) stacking , as shown in fig . [ fig_schem ] . a unit cell contains @xmath6 and @xmath7 atoms on the layer @xmath8 . for convenience we divide carbon atoms into two groups as @xmath9 the atoms of group i are arranged along vertical columns normal to the layer plane , while those in group ii are above or below the center of hexagons in the neighboring layers . the lattice constant within a layer is given by @xmath10 nm and the distance between adjacent layers @xmath11 nm . the system can be described by a * k*@xmath12*p * hamiltonian closely related to a three - dimensional ( 3d ) graphite model.@xcite the low energy spectrum is given by the states in the vicinity of @xmath1 and @xmath2 points in the brillouin zone . let @xmath13 and @xmath14 be the bloch functions at the @xmath1 point , corresponding to the @xmath15 and @xmath16 sublattices , respectively , of layer @xmath17 . for monolayer graphene , the hamiltonian around @xmath1 point for the basis @xmath18 , @xmath19 is written as @xcite @xmath20 where @xmath21 and @xmath22 with @xmath23 being the nearest - neighbor coupling in a single layer . we cite the experimental estimation @xmath24 ev.@xcite    for the inter - layer coupling , we include parameters @xmath25 , and @xmath26 , following the hamiltonian previously derived for a bilayer graphene.@xcite here @xmath25 represents the coupling between vertically neighboring atoms in group i ( @xmath27 ) , and @xmath26 between group ii atoms on neighboring layers ( @xmath28 ) , which are estimated to @xmath29 ev,@xcite and @xmath30 ev.@xcite if we look at the interaction between layers 1 and 2 , the matrix element @xmath31 , corresponding to the vertical bond , becomes @xmath25 not accompanying in - plane bloch number . the matrix element @xmath32 is written as @xmath33 with @xmath34 similar to the intra - layer term @xmath35 , as the in - plane vector components from @xmath36 to @xmath37 are identical to those from @xmath38 to @xmath36 . accordingly , if the basis is taken as @xmath39 ; @xmath40 ; @xmath41 ; @xmath42 , the hamiltonian for the multilayer graphene around the @xmath1 point becomes @xmath43 with @xmath44 the effective hamiltonian for @xmath2 is obtained by exchanging @xmath45 and @xmath46 and replacing @xmath25 with @xmath47 . the derivation of the effective hamiltonian based on a tight - binding model is presented in appendix [ sec_app_a ] . we show in the following that the hamiltonian matrix ( [ eq_h ] ) can be block - diagonalized into smaller matrices by choosing an appropriate basis independent of @xmath48 . first , we arrange the basis in the order of group i and then group ii , i.e. , @xmath19 , @xmath49 , @xmath50 , @xmath41 ; @xmath18 , @xmath51 , @xmath52 , @xmath41 . then , eq .  ( [ eq_h ] ) becomes @xmath53 with @xmath54 being @xmath55 matrices defined as @xmath56 where the upper and lower signs correspond to odd and even @xmath5 , respectively .    if we set @xmath57 ( the @xmath1 point ) , @xmath58 and @xmath59 vanish . remaining @xmath60 is equivalent to the hamiltonian of a one - dimensional tight - binding chain with the nearest - neighbor coupling @xmath25 , giving a set of eigenenergies @xmath61 with @xmath62 here , @xmath63 is an odd integer when the layer number @xmath5 is even , while @xmath63 is even when @xmath5 is odd , and therefore @xmath64 is allowed only for odd @xmath5 .    the corresponding wave function is explicitly written as @xmath65 , \label{eq_psi}\ ] ] where @xmath66 represents the amplitudes at @xmath19 , @xmath49 , @xmath50 , @xmath41 and satisfies @xmath67 we have a relation between the wave functions @xmath68 and @xmath69 as @xmath70 now we construct the basis by assigning @xmath71 to the atoms of group i and ii as @xmath72 and attempt to rewrite the hamiltonian ( [ eq_h ] ) . the matrix elements within group i come from @xmath60 and become diagonal as is obvious from the definition , @xmath73 off - diagonal elements between @xmath74 and @xmath75 are written from @xmath58 as @xmath76 in the second equality we used relation ( [ eq_bipartite ] ) and orthogonality ( [ eq_orthogonality ] ) . lastly , the matrix elements within group ii are obtained from @xmath59 as @xmath77 \nonumber \\ + \!\ ! & \!\!\!\ ! & \gamma ' i k_y \sum_{j=1}^{n-1 } ( -1)^j   \left [ \psi_{m'}^*(j+1 ) \psi_{m}(j ) - \psi_{m'}^*(j ) \psi_{m}(j+1 )   \right ] \nonumber\\ = & \!\!\!\ ! & \gamma ' \lambda_{n , m } ( k_x \delta_{m , m ' } + i k_y \delta _ { m ,- m ' } ) . \ ] ]    the hamiltonian is thus closed in the subspace @xmath78 for each @xmath79 . particularly , @xmath64 is special in that the subspace is spanned with only two bases @xmath80 , while this is absent in even - layer graphenes . the sub - matrix is written as @xmath81 which is independent of @xmath25 and @xmath26 , and equivalent to the hamiltonian of the monolayer graphene .    for @xmath82 , we rearrange the basis as @xmath83 where we take @xmath84 without loss of generality . we then obtain @xmath85 with @xmath86 . this is equivalent to the hamiltonian of a bilayer graphene except that @xmath25 and @xmath87 @xmath88 are multiplied by @xmath89 . thus the hamiltonian of odd - layered graphene is composed of one monolayer - type and @xmath90 bilayer - type subbands while that of even - layered graphene is composed of @xmath91 bilayers but no monolayer . the similar idea was previously proposed for trilayer graphene without @xmath26 , where it was shown that the energy spectrum becomes a superposition of that for a monolayer and for a bilayer.@xcite here we have extended this argument to decomposition of the hamiltonian matrix , and to systems with arbitrary number of layers including the trigonal warping . we also note that @xmath48-independence of the basis becomes important in the following sections , since this enables us to write the magnetization as a sum over contributions from sub - hamiltonians , which are independently calculated . many other parameters were introduced for the description of the band structure of bulk graphite.@xcite the parameter @xmath92 couples group i and ii atoms sitting on the neighboring layers , such as @xmath93 or @xmath94 . this parameter does not change the qualitative feature of the low - energy spectrum and therefore is not important.@xcite parameters @xmath95 and @xmath96 represent vertical hoppings between the second - nearest neighboring layers for group ii and i atoms , respectively . further , @xmath97 is an energy difference between the group i and ii atoms due to difference in the chemical environment . inclusion of these parameters @xmath95 , @xmath96 , and @xmath97 causes opening up of small energy gaps between the conduction and the valence bands . however , these gaps do not play important roles in the magnetization as will be discussed in the following .    in 3d limit , @xmath98 , the eigenstate becomes a superposition of opposite traveling waves with @xmath99 along the stacking direction . the relation between the index @xmath63 and @xmath100 is obtained by comparing the eigenenergy of @xmath60 , eq . ( [ eq_1d ] ) , to that of the 3d limit , @xmath101 , as @xmath102    the band structure of the hamiltonian ( [ eq_hm ] ) can be obtained by replacing @xmath25 by @xmath103 and @xmath26 by @xmath104 in that of the bilayer.@xcite we plot in fig . [ fig_band0 ] the dispersion for @xmath105 , which has the maximum trigonal warping . the middle two subbands stick together at @xmath106 while the remaining two bands appear only in the energy range @xmath107 . if we neglect @xmath26 , the effective hamiltonian for @xmath108 becomes @xmath109 which works for the reduced basis @xmath110 giving a rotationally symmetric dispersion with the effective mass @xmath111    ) with @xmath105 and @xmath112 . @xmath113 is shown as horizontal dotted lines . right panel shows zoom out of the left . ( bottom ) 3d plot of the lower second band around the band touching point . four fermi points at @xmath4 indicated by dots . , width=302 ]    the term proportional to @xmath26 is responsible for the trigonal warping effect , which is most remarkable around the band sticking point @xmath114 . let us define @xmath115 in the energy range @xmath116 , the fermi line splits into four separated pockets , one center part and three leg parts located trigonally , which shrink into four fermi points linearly with @xmath117 . we note that @xmath118 is proportional to @xmath119 and thus very sensitive to @xmath89 , while the energy of the higher - band bottom , @xmath103 , behaves linear to @xmath89 . the maximum of @xmath89 approaches 2 as the layer number increases , so that @xmath118 becomes as large as @xmath120 mev . figure [ fig_band ] shows the band structures around the @xmath1 point along the @xmath121 axis in the multilayer graphenes with @xmath122 , 3 , 4 , and 5 and @xmath123 . the lists of @xmath124 are given as @xmath125    while we have included @xmath23 , @xmath25 , and @xmath26 in our graphene model , the extra parameter neglected here may make some changes in the electronic structure . the energy band of a few - layered graphene has been calculated in the density functional calculation @xcite and the tight - binding model.@xcite those results differ from ours mainly in that the band centers relatively shift depending on @xmath63 , and that a narrow gap opens where the conduction and valence bands ( within a single @xmath63 ) touch , and where different bands ( with different @xmath63 s ) cross . gaps are attributed to effects of couplings such as @xmath95 , @xmath96 , and @xmath97 , which are mentioned above . in terms of the effective mass hamiltonian ( [ eq_h ] ) , those parameters appear as matrix elements without being multiplied by the wave number @xmath121 and @xmath126 , since they are associated with a hopping along the @xmath127 axis or a diagonal element . thus , they do not vanish at @xmath128 ( @xmath1 or @xmath2 ) and lift the degeneracy to open a gap . apart from the gap opening , the main feature of the trigonal warping is well described in the present model . it should also be mentioned that an energy gap is induced by an electric field perpendicular to the layer stacking direction,@xcite where the electrostatic potential appears as matrix elements independent of @xmath121 and @xmath126 as well .     with @xmath123 , around the @xmath1 point ( taken as origin ) along the @xmath121 axis . right panel shows a zoom - out of the left . numbers assigned to curves indicate @xmath63 . , width=302 ] for the magnetic susceptibility , we use the general expression based on the linear response theory,@xcite @xmath129 with @xmath130 where @xmath131 is the valley degeneracy , @xmath132 is the spin degeneracy , and @xmath133 is the system size . we defined here @xmath134 , @xmath135 , @xmath136 , and @xmath137^{-1}$ ] with the chemical potential @xmath138 and the temperature @xmath139 . the formula valid also for the hamiltonian ( [ eq : bilayer_low_energy ] ) is discussed in appendix [ sec_app_b ] . by integration by parts in eq . ( [ eq_chi ] ) , we have @xmath140 showing that the susceptibility at non - zero temperature is written in terms of that at zero temperature . the integration of @xmath141 over @xmath138 is independent of @xmath139 . we include the impurity scattering effects by introducing a self - energy @xmath142 in the green s function , i.e. , @xmath143 in ( [ eq_chi ] ) is replaced by @xmath144 . here we simply assume the scattering rate @xmath145 to be independent of energy . using the decomposition of the hamiltonian , the magnetization of the @xmath5-layered graphene can be written as a summation over each sub - hamiltonian . the contribution from @xmath64 is exactly equivalent to the susceptibility of a monolayer graphene,@xcite which becomes at zero temperature and in the clean limit , @xmath146 thus , the odd - layer graphene always has a large diamagnetic peak at zero energy . the delta - function dependence of @xmath147 agrees with the general property of the susceptibility in systems described by the @xmath148-linear hamiltonian , as discussed in sec . [ sec_disc ] .    in the presence of disorder , the delta - function is broadened into a lorentzian with width @xmath149 and the same area,@xcite i.e. @xmath150 within the present model assuming a constant @xmath149 . the shape of the peak itself depends on the model disorder and we may have some different manner of broadening in a more realistic treatment . in fact , in the monolayer graphene it was shown in a self - consistent born approximation@xcite that @xmath141 has a much sharper peak at @xmath151 than the lorentzian and also a large tail proportional to @xmath152 for @xmath153.@xcite in multi - layer cases effects of disorder are more complicated because of the presence of other bands . this problem is out of the scope of this work . the susceptibility of a bilayer graphene described by the hamiltonian ( [ eq_hm ] ) was analytically calculated for the case of @xmath154.@xcite the expression for @xmath155 and @xmath156 is given by @xmath157 with @xmath86 , where @xmath158 is a step function defined by @xmath159 the susceptibility diverges logarithmically toward @xmath0 , becomes slightly positive for @xmath160 , and vanishes for @xmath161 where the higher subband enters . in the presence of disorder , the logarithmic peak is broadened approximately as @xmath162 . the integration of @xmath141 in eq . ( [ eq_chi_bi ] ) over the fermi energy becomes @xmath163 independent of @xmath25 , which is exactly twice as large as that of the monolayer graphene ( [ eq_chi_mono ] ) . this arises due to the fact that the integral of @xmath141 over the fermi energy is determined only by terms of the hamiltonian matrix , proportional to @xmath121 or @xmath126 , and is independent of terms independent of @xmath121 and @xmath126 . a proof of this important property is presented in sec . [ sec_disc ] .    if we include the extra band parameter @xmath26 , the low - energy structure of the susceptibility ( [ eq_chi_bi ] ) drastically changes due to the fine structure around the band touching point . to demonstrate this , we numerically calculate @xmath141 for the hamiltonian ( [ eq_hm ] ) in the case of the maximum trigonal warping , @xmath164 . figure [ fig_chi ] shows the susceptibility as a function of @xmath165 with several values of @xmath149 . we take @xmath166 , where the fermi line splitting occurs in lower than @xmath167 . for reference we also plot the result without the trigonal warping , @xmath168 , as a dashed curve .    ) in the multilayer graphene with @xmath169 , @xmath123 and several disorder strengths @xmath149 . energy is scaled in units of @xmath167 . dashed curves show plots for @xmath170 . inset at the top is a zoom out of the top panel ( @xmath171 ) with units of energy @xmath25 . , width=302 ]    when we go down from high energy in the top panel ( the smallest @xmath149 ) , the susceptibility gradually deviates downward from the logarithmic dependence of @xmath168 , and takes a sharp dip at @xmath172 . remarkably we have a strong peak centered on @xmath4 , which is regarded as the effect of the linear dispersions around zero energy . the integral of @xmath141 over the fermi energy is almost constant @xmath173 as discussed in sec . [ sec_disc ] , showing that the reduction in higher energies compensates the zero energy peak . as @xmath149 becomes larger , the peak begins to cancel with the reduction in high energy and the effect of @xmath26 eventually disappears when @xmath174 . in contrast , the peak associated with the monolayer band @xmath64 becomes broad in @xmath149 but never vanishes , as shown in eq.([eq : monolayer_with_broadening ] ) . figure [ fig_chi_multi ] shows @xmath175 of graphenes with layer number from @xmath176 to 5 with several disorder strengths @xmath149 . for @xmath177 , insets show the contributions from each of bilayer - type bands . the result of odd @xmath5 always contains a monolayer - like component , which is exactly the same as @xmath176 and thus omitted in the inset . we can see that odd - layered graphenes exhibit a particularly large peak , which mainly comes from the monolayer - type band . a bilayer - like component contains a central peak due to the trigonal warping and a logarithmic tail in high energies , in accordance with fig . [ fig_chi ] . the layer - number dependence of the susceptibility in multilayer graphene has been studied for the graphite intercalation compounds.@xcite this system can be viewed as independent multilayer graphenes bound by the intercalant layers , but the intercalants give a strong electrostatic potential along the stacking direction , leading to the charge redistribution among different layers.@xcite as a result the band structure and the magnetization are considerably different from our system with a uniform electrostatic potential in the vertical direction .    in isolated multilayer graphenes realized in recent experiments , we may have some potential difference among layers depending on the experimental environment , and this can also be tuned by the external electric field as mentioned . in sec .  [ sec_disc ] , we will show that , as long as the potential is not too strong to alter the entire band structure , this does not change the qualitative feature of the magnetization . the zero energy peak in the bilayer - type subband originates in dirac - like dispersions appearing around four fermi points . using the known results in a bilayer,@xcite we can show that the sequence of the landau levels in the center pocket approximately becomes @xmath178 with @xmath179 , and those in the three leg parts @xmath180 , where @xmath181 is the magnetic length . since the susceptibility is determined solely by landau level energies , we compare this to the monolayer s sequence @xmath182 and obtain @xmath141 from each pocket by substituting @xmath183 in eq . ( [ eq_chi_mono ] ) . we end up with @xmath184 except for a constant coming from the integral over the lower energy states . the zero - energy peak in fig . [ fig_chi ] fits well to the lorentzian with width @xmath149 and the area of the delta - function ( [ eq_chi_fine ] ) , as long as @xmath185 . the factor attached to @xmath147 becomes as large as 0.4 when @xmath164 and @xmath123 , and therefore the singularity is not too small compared with that of the monolayer . for @xmath5-layered graphene , a simple relation @xmath186 leads to the summation of ( [ eq_chi_fine ] ) over all the bilayer - type subbands , @xmath187 in fig . [ fig_chi_multi ] the peak height becomes a little larger than this estimation due to mixing with the logarithmic tail . the delta - function dependence of @xmath141 in monolayer graphene is a characteristic property common to general @xmath148-linear hamiltonian . this can be shown using the scaling argument . we consider a hamiltonian @xmath188 which contains only terms linear in @xmath121 and @xmath126 . we change the energy and wave number scales by an arbitrary factor @xmath189 as @xmath190 then the hamiltonian becomes formally identical under this transformation , since the coefficients of @xmath148-linear terms in the hamiltonian remain unchanged . going back to the definition of @xmath141 in ( [ eq_chi ] ) and ( [ eq_f ] ) , @xmath191 is scaled as @xmath192 the function @xmath193 should depend only on the coefficients of @xmath148-linear terms and natural constants , and thus is invariant under the scale transformation , namely we have @xmath194 . with ( [ eq_f_scale ] ) , we come up with a equation , @xmath195 which is satisfied solely by @xmath196    a constant @xmath15 is related to the integral of the susceptibility @xmath175 over the fermi energy @xmath165 . from ( [ eq_chi ] ) , we generally have @xmath197 where the integral path @xmath198 is a circle with an infinite radius with anti - clockwise direction . in the present system , ( [ eq_f_explicit ] ) immediately gives the integral as @xmath199 . this is an integral of the real function @xmath200 and thus is real . substituting eq . ( [ eq_f_explicit ] ) with real @xmath15 in eq . ( [ eq_chi ] ) , we finally obtain the explicit form of the zero - temperature susceptibility as @xmath201    as discussed in sec . [ sec_form ] , the band structure in more realistic models has an energy gap around zero energy due to extra band parameters neglected in the present model . it was also mentioned that the external electric field along stacking direction opens an energy gap . one might think that the gap would strongly reduce the large diamagnetism at the band touching point . however , we can show within the effective mass approximation that the integral of susceptibility over @xmath165 is independent of any kind of matrix elements without @xmath121 and @xmath126 , which are responsible for gap opening . this is obvious from the general expression ( [ eq_chi_oint ] ) ; even if the hamiltonian contains @xmath148-independent terms in addition to @xmath148-linear terms , they can be safely neglected in the integral as they are infinitesimal compared to @xmath202 on the path @xmath198 . in the effective - mass model of the multilayer graphene we immediately conclude that the integral is independent of @xmath25 , @xmath95 , @xmath96 , @xmath97 , and any other parameters independent of the wave vector . thus we expect that the large diamagnetic peak is still visible even when a gap opens , while it may get broadened in energy by the gap width . the diamagnetism in narrow gap systems is known in bismuth@xcite and recently studied for the gapped dirac fermion.@xcite any further discussion requires a direct computation of the magnetization including extra parameters , but we leave this for the future study . the integral of the susceptibility can be calculated by the hamiltonian with the @xmath148-independent terms dropped , and thus depends only on the band parameters associated with @xmath148-linear terms . in our model ( [ eq_h ] ) , the value is mainly determined by the dominant parameter @xmath23 , while @xmath26 gives a correction at most of the order of @xmath203 . the correction must be the second order in @xmath26 because we can change @xmath26 to @xmath204 in the hamiltonian with a unitary transformation multiplying the base on layer @xmath17 by @xmath205 . as a result , the integral of @xmath175 for the bilayer - type hamiltonian becomes almost twice as large as the monolayer s , and the summation over all the subsystems in @xmath5-layered graphene becomes approximately @xmath5 times as large as the monolayer s . it is instructive to derive the susceptibility starting from the landau - level energies . in the monolayer graphene , the thermodynamic function @xmath206 is given by @xmath207 \big\ } , \end{aligned}\ ] ] where @xmath208 , @xmath209 with @xmath210 , and @xmath211 a cutoff function which gradually decays to zero for @xmath212 with cutoff energy @xmath213 . we can rewrite this as @xmath214 where @xmath215 , \ ] ] with @xmath216 .    expanding the integral @xmath217 with respect to @xmath218 , we immediately have @xmath219 \nonumber \\ = & \!\!\!\!\!\ ! & \int_0^\infty h(x ) d x - { 1\over 12 } h^2 \big [ h'(0 ) \!+\ ! { 1\over 2 } h'(\infty ) \big ] , \\end{aligned}\ ] ] up to the second order in @xmath218 or in @xmath16 . then , we have @xmath220 where @xmath221 is the thermodynamic function in the absence of a magnetic field and @xmath222 ^ 2 } \nonumber\\ & = & { g_{\rm v } g_{\rm s } \gamma^2 \over 12 \pi l^4 } \int_{-\infty}^{\infty } \big ( \!-\ ! { \partial f(\varepsilon ) \over \partial \varepsilon } \big ) \delta(\varepsilon ) d \varepsilon . \label{eq : thermodynamic_potential:_monolayer_graphene}\end{aligned}\ ] ] applying the relation @xmath223 , we obtain @xmath147 given by ( [ eq_chi_mono ] ) at zero temperature . we should note that the thermodynamic function in the absence of a magnetic field is given by @xmath224 for contributions of states with @xmath225 , we can expand the above with respect to @xmath226 and have to the lowest order in the field strength @xmath16 @xmath227 \nonumber\\ & \ ! \approx \ ! & { g_{\rm v } g_{\rm s } \over 2\pi l^2 } { ( \hbar\omega_b)^2 \over 48 } \!\ ! \lim_{\delta\rightarrow+0 } \big ( \int_\delta^\infty [ \varepsilon^{-2 } f(\varepsilon ) \!-\ ! \varepsilon^{-1 } f'(\varepsilon ) ] d \varepsilon \nonumber \\ & & \qquad\qquad\qquad\quad   - \int_{-\infty}^{-\delta } [ \varepsilon^{-2 } f(\varepsilon ) \!-\ ! \varepsilon^{-1 } f'(\varepsilon ) ] d \varepsilon \big ) \nonumber \\ & \ ! = \ ! & - { g_{\rm v } g_{\rm s } \over 2\pi l^2 } { 1\over 24 } ( \hbar\omega_b)^2 \int_{-\infty}^\infty \big ( \!-\ ! { \partial f(\varepsilon ) \over \partial \varepsilon } \big ) \delta(\varepsilon ) d \varepsilon .\end{aligned}\ ] ] this gives a `` paramagnetic '' susceptibility . for @xmath228 , on the other hand , the change in the thermodynamic potential is calculated as @xmath229 the sum of these two contribution is the same as eq . ( [ eq : thermodynamic_potential:_monolayer_graphene ] ) , as is expected .    in the bilayer graphene , the landau level with @xmath154 in the region @xmath230 can be calculated from the hamiltonian ( [ eq : bilayer_low_energy ] ) as @xcite @xmath231 where @xmath232 with @xmath233 defined in eq . ( [ eq_mstar ] ) , @xmath234 , and @xmath235 . we have doubly degenerate levels at zero energy ( @xmath236 , @xmath234 ) , while the spacing gradually becomes constant as @xmath237 goes higher . in a similar but more complicated manner , the susceptibility is calculated as @xmath238 which correctly describes the logarithmic divergence around zero energy in the rigorous expression ( [ eq_chi_bi ] ) , as expected . a constant term independent of energy is missing in eq . ( [ eq_chi_bi_alt ] ) since this depends on all the low - energy bands which are neglected in this calculation . we can understand the logarithmic dependence intuitively by looking into the landau - level sequence . the landau level energy can be expanded for large @xmath237 as @xmath239,\ ] ] where the first term gives the constant interval , and the second gives a shift toward zero energy , which is rewritten as @xmath240 . for @xmath241 , for example , the change in the total energy due to the energy shift is calculated as @xmath242 giving the @xmath243 dependence of the susceptibility .    for the hamiltonian ( [ eq : bilayer_low_energy ] ) containing terms proportional to @xmath244 , the susceptibility formula ( [ eq_chi ] ) with ( [ eq_f ] ) is no longer valid , since this was originally derived for systems in which @xmath245 commutes with @xmath246 and @xmath247 . the modified formula should be @xmath248 this is derived in appendix [ sec_app_b ] . a scaling argument similar to the case of the monolayer graphene then gives @xmath249 . this again leads to the logarithmic dependence of @xmath141 on the fermi energy , which coincides with ( [ eq_chi_bi_alt ] ) apart from a constant . the experimental measurements of the magnetization of two - dimensional electron systems were performed on the semiconductor heterostructures , by using the superconducting quantum interference device ( squid)@xcite or using the torque magnetometer.@xcite we expect that the detection of the graphene magnetism is also feasible with those techniques .    we have studied the orbital magnetism of multilayer graphene with the bernal stacking in the effective mass approximation . we have demonstrated that the hamiltonian and thus the susceptibility can be decomposed into those equivalent to the monolayer or bilayer bands . the monolayer - like band exists only in odd - layered graphenes and gives a strong diamagnetic peak at @xmath250 . the bilayer - like bands always exist and present a strong diamagnetism in the vicinity of zero energy , unless the fine band structure caused by @xmath26 is destroyed by the disorder . this work has been supported in part by the 21st century coe program at tokyo tech nanometer - scale quantum physics  and by grants - in - aid for scientific research from the ministry of education , culture , sports , science and technology , japan . we derive in the following the effective mass equation eq . ( [ eq_h ] ) describing states in the vicinity of @xmath1 point in a multilayer graphene , by starting from the one - orbital tight - binding model . the following is nothing but a straightforward extension of the monolayer case.@xcite in a tight - binding model , the wave function is written as @xmath251,\end{aligned}\ ] ] where @xmath252 is the layer index , @xmath253 is the wave function of the @xmath254 orbital of a carbon atom located at the origin , as a function of three - dimensional position @xmath255 . @xmath256 is the three - dimensional position of the site @xmath257 , and @xmath258 is a two - dimensional component of @xmath259 , parallel to the layer .    in the model including hopping parameters @xmath23 , @xmath25 , and @xmath26 defined in sec . [ sec_form ] , schrdinger s equation can be written as follows : for odd @xmath17 , @xmath260 , \\   \vare \psi_{b_{j}}({\mbox{\boldmath$\rho$}}_{b_{j } } ) & = & -\gamma_0 \sum_{l=1}^3    \psi_{a_{j}}({\mbox{\boldmath$\rho$}}_{b_{j } } + { \mbox{\boldmath$\tau$}}_l ) \nonumber\\ & & \hspace{0 mm } + \gamma_1 \left [ \psi_{a_{j+1}}({\mbox{\boldmath$\rho$}}_{b_{j } } ) + \psi_{a_{j-1}}({\mbox{\boldmath$\rho$}}_{b_{j } } ) \right ] . \quad \label{eq_sch1}\end{aligned}\ ] ] for even @xmath17 , @xmath261 , \quad \\   \vare \psi_{b_{j}}({\mbox{\boldmath$\rho$}}_{b_{j } } ) & = & -\gamma_0 \sum_{l=1}^3    \psi_{a_{j}}({\mbox{\boldmath$\rho$}}_{b_{j } } + { \mbox{\boldmath$\tau$}}_l ) \nonumber\\ & & \hspace{-20 mm } + \gamma_3 \sum_{l=1}^3    \left [ \psi_{a_{j+1}}({\mbox{\boldmath$\rho$}}_{b_{j } } - { \mbox{\boldmath$\tau$}}_l ) + \psi_{a_{j-1}}({\mbox{\boldmath$\rho$}}_{b_{j } } - { \mbox{\boldmath$\tau$}}_l ) \right ] . \label{eq_sch2}\end{aligned}\ ] ] here we introduced the vectors from b site to the nearest neighboring a sites as @xmath262 , @xmath263 , and @xmath264 , and we set @xmath265 .    the states around @xmath1 point can be expressed in terms of the slowly - varying envelope functions @xmath266 as @xmath267 where @xmath268 , and @xmath269 , @xmath270 are phase factors defined by @xmath271 with @xmath272 . when @xmath273 is much smaller than the length scale of the envelope functions , we have @xmath274 with @xmath275 or @xmath7 .    by substituting eq . ( [ eq_env ] ) with ( [ eq_approx ] ) , into schrdinger s equations ( [ eq_sch1 ] ) and ( [ eq_sch2 ] ) , we have for odd @xmath17 , @xmath276 \nonumber\\   \vare f_{b_j}({\mbox{\boldmath$\rho$ } } ) & = &     \gamma k_+   f_{a_j}({\mbox{\boldmath$\rho$ } } ) + \gamma_1   \left [ f_{a_{j-1}}({\mbox{\boldmath$\rho$ } } ) + f_{a_{j+1}}({\mbox{\boldmath$\rho$ } } ) \right ] , \nonumber\\\end{aligned}\ ] ] and for even @xmath17 , @xmath277 , \nonumber\\   \vare f_{b_j}({\mbox{\boldmath$\rho$ } } ) & = &     \gamma k_+   f_{a_j}({\mbox{\boldmath$\rho$ } } ) + \gamma ' k_-    \left [ f_{a_{j-1}}({\mbox{\boldmath$\rho$ } } ) + f_{a_{j+1}}({\mbox{\boldmath$\rho$ } } ) \right],\nonumber\\\end{aligned}\ ] ] where @xmath21 with @xmath278 , @xmath279 . we also used an identity @xmath280 if we rewrite this set of equations into the matrix form for a vector @xmath281 , we finally obtain the hamiltonian matrix ( [ eq_h ] ) .    the effective hamiltonian for another valley @xmath282 can be derived in a parallel way , while @xmath283 and @xmath284 are exchanged in eq . ( [ eq_c ] ) . the susceptibility formula ( [ eq_chi ] ) with ( [ eq_f ] ) has been derived for the luttinger - kohn representation of the bloch function@xcite and therefore in systems described by the hamiltonian consisting of the free electron kinetic energy @xmath285 and terms linear in @xmath286.@xcite we derive here the susceptibility formula which is valid in the general hamiltonian @xmath287 which includes @xmath148-square terms in off - diagonal matrix elements as well as @xmath148-linear terms . we shall confine ourselves to the case without electron - electron interaction for simplicity . consider the system described by the schrdinger equation @xmath288 with @xmath289 and @xmath290 being the vector potential . the thermodynamic function @xmath206 is given by @xmath291 \big\ }   \nonumber\\   & = & - k_{\rm b } t \ , g_{\rm s } \frac{1}{v }   \int d\varepsilon \ , \big(\!-\!{1\over\pi}\big ) { \rm i m } \ , { \rm tr } \ ,   \frac{1}{\vare -{\cal h } + i0 } \nonumber\\ & & \qquad\qquad \times\ln \big\ { 1 \!+\ ! \exp[\beta(\mu \!-\ ! \varepsilon ) ] \big\ } ,   \label{eq_omega}\end{aligned}\ ] ] where @xmath292 is the spin degeneracy and @xmath293 is the system volume . we consider an isotropic system and assume the vector potential @xmath294 , with @xmath295 where we are going to take the long wavelength limit @xmath296 , for which the field causes the response the same as that due to a spatially uniform magnetic field . in the presence of this vector potential , the hamiltonian changes from @xmath297 to @xmath298 , with @xmath299 , where @xmath300 . the hamiltonian can be expanded as @xmath301 where @xmath302 , @xmath303 , and @xmath304 . note that in general @xmath148-square hamiltonian @xmath305 does not commute with @xmath246 but does with @xmath306 . expanding the hamiltonian up to the second order in the strength of the magnetic field @xmath16 , we have @xmath307 where @xmath308 is the magnetic length , @xmath309 , @xmath310 , and we assumed that the system is translational invariant ( after the configuration average in the presence of impurities ) . we then expand @xmath311 up to the second order in @xmath312 , to have @xmath313 with @xmath314 . this immediately gives the change of the thermodynamic potential @xmath315 with eq . ( [ eq_omega ] ) . the susceptibility @xmath141 is obtained by a relation @xmath316 as @xmath317 this gives eq . ( [ eq_chi ] ) with ( [ eq_f_gen ] ) for the multilayer graphene . d. p. divincenzo and e. j. mele , phys . b * 29 * , 1685 ( 1984 ) . j. gonzalez , f. guinea , and m. a. h. vozmediano , nucl . b * 424 * , 595 ( 1994 ) . j. gonzalez , f. guinea , and m. a. h. vozmediano , phys . lett .  * 77 * , 3589 ( 1996 ) . n. h. shon and t. ando , j. phys . * 67 * , 2421 ( 1998 ) . j. gonzalez , f. guinea , and m. a. h. vozmediano , phys . b * 63 * , 134421 ( 2001 ) . y. zheng and t. ando , phys . b * 65 * , 245420 ( 2002 ) . t. ando , y. zheng , and h. suzuura , j. phys . jpn .  * 71 * , 1318 ( 2002 ) . h. suzuura and t. ando , phys . lett .  * 89 * , 266603 ( 2002 ) . t. ando , j. phys . * 74 * , 777 ( 2005 ) . k. s. novoselov , e. mccann , s. v. morozov , v. i. falko , m. i. katsnelson , u. zeitler , d. jiang , f. schedin , and a. k. geim , nat . phys . * 2 * , 177 ( 2006 ) . t. ohta , a. bostwick , t. seyller , k. horn , and e. rotenberg , science * 313 * , 951 ( 2006 ) ; t. ohta , a. bostwick , j. l. mcchesney , t. seyller , k. horn , and e. rotenberg , phys . lett . * 98 * , 206802 ( 2007 ) . e. mccann and v. i. falko , phys . lett . * 96 * , 086805 ( 2006 ) .
we present a theoretical study on the orbital magnetism in multilayer graphenes within the effective mass approximation . 
 the hamiltonian and thus susceptibility can be decomposed into contributions from sub - systems equivalent to monolayer or bilayer graphene . 
 the monolayer - type subband exists only in odd layers and exhibits a delta - function susceptibility at @xmath0 . 
 the bilayer - type subband appearing in every layer number gives a singular structure in the vicinity of @xmath0 due to the trigonal warping as well as a logarithmic tail away from @xmath0 . 
 the integral of the susceptibility over energy is approximately given only by the layer number .