document: recently , topological properties of time - reversal - invariant band insulators in two and three dimensions have been extensively studied@xcite . a class of insulators preserving the time reversal symmetry is called topological insulators characterized by non - trivial topological invariants@xcite.the topological insulators 
   have been intensively studied because of the existence and potential applications of robust surface metallic states .    both in two and three dimensions , the topological phases are typically realized in the systems with strong spin - orbit interaction@xcite . 
   all the known topological insulators contain heavy or rare metal elements , such as bismuth or iridium , which poses constraints on the search for topological materials . 
   irrespective of constitutents , ubiquitous mutual coulomb repulsions among electrons have been proposed to generate effective spin - orbit couplings  @xcite . 
   it has been proposed that an extended hubbard model on the honeycomb lattice can generate an effective spin - orbit interaction from a spontaneous symmetry breaking at the hartree - fock mean - field level leading to a topologically non - trivial phase@xcite . since the honeycomb - lattice system , which is dirac semimetals in the non - interacting limit , becomes a topologically nontrivial insulator driven by the coulomb interaction , this phase 
   is often called a topological mott insulator ( tmi ) . 
   this phenomenon is quite unusual not only because an emergent spin - orbit interaction appears from the electronic mutual coulomb interaction , but also it shows an unconventional quantum criticality that depends on the electron band dispersion near the fermi point@xcite . 
   however , this proposed topological phase by utilizing the ubiquitous coulomb repulsions has not been achieved in real materials even though the tmi is proposed not only in various solids  @xcite but also in cold atoms loaded in optical lattices  @xcite . even in simple theoretical models such as extended hubbard models , 
   it is not clear whether the tmis become stable against competitions with other orders and quantum fluctuations . 
   reliable examination of stable topological mott orders in the extended hubbard model is hampered by competing symmetry breakings such as cdws . 
   couplings driving the topological mott transitions are also relevant to formations of a cdw , which has not been satisfactorily discussed in the previous studies . since the emergence of the tmi in the honeycomb lattice requires the coulomb repulsion between the next nearest neighbor sites , the long - period cdw instability must be considered on equal footing , which is not captured in the small - unit - cell mean - field ansatz employed in the previous studies . 
   examining charge fluctuations with finite momentum over entire brillouin zones is an alternative way to clarify the competitions among tmis and cdws , as studied by employing functional renormalization group methods  @xcite . 
   however , first order thermal or quantum phase transitions not characterized by diverging order - parameter fluctuations are hardly captured by such theoretical methods . 
   the most plausible symmetry breking competing with tmis indeed occurs as a first order quantum phase transition as discussed later . 
   the quantum many - body fluctuations beyond the mean - field approximation severely affects the stability of the tmi . 
   the stability of the tmi and estimation of the critical value of interaction on the honeycomb lattice has mainly been considered by mean - field calculations which can not treat the correlation effect satisfactorily . 
   however , there exists a reliable limit where the tmi becomes stable : for infinitesimally small relevant coulomb repulsions , the quadratic band crossing with vanishing fermi velocities cause the leading instability toward the tmi , as extensively examined by using perturbative renormalization group methods@xcite . 
   however , examining the instabilities toward the tmi in dirac semimetals requires elaborate theoretical treatments .    in this study , for clarification of the competitions among tmis and other symmetry breakings 
   , we first examine the long - period cdw at the level of mean - field approximation that turns out to be much more stable compared to that of short period . indeed , this cdw severly competes the tmi on the honeycomb lattice . 
   the tmi on the honeycomb lattice studied in the literatures is consequently taken over by the cdw .    we , however , found a prescription to stabilize the tmis on the honeycomb lattice : by reducing the fermi velocity of the dirac cones , the tmi tends to be stabilized . 
   we examine the realization of the tmis in the extended hubbard model on the honeycomb lattice by controlling the fermi velocity and employing a variational monte carlo method@xcite with many variational parameters@xcite , multi - variable variational monte carlo ( mvmc)@xcite , together with the mean - field approximation . 
   finally , we found that , by suppressing the fermi velocity to a tenth of that of the original honeycomb lattice , the tmi emerges in an extended parameter region as a spontaneous symmetry breaking even when we take many - body and quantum fluctuations into account . 
   this paper is organized as follows . in section 
   [ sec : model and method ] , we introduce an extended hubbard model and explain the order parameter of tmi . 
   we also introduce the mvmc method . in section [ sec : stability ] , we first show how the long - range cdw becomes stable over the tmi phase in standard honeycomb lattice models . 
   then we pursue the stabilization of tmi by modulating fermi velocity at the dirac cone at the mean - field level . 
   finally we study by the mvmc method the effect of on - site coulomb interaction which was expected to unchange the stability of the tmi phase at the level of mean - field approximation . 
   section [ sec : dis ] is devoted to proposal for realization of tmis in real materials such as twisted bilayer graphene . 
   in this section , we study ground states of an extended hubbard model on the honeycomb lattice at half filling defined by @xmath0 where the single particle parts of @xmath1 are defined as @xmath2 and @xmath3 is the spin - orbit interaction . here @xmath4 is a creation ( annihilation ) operator for a @xmath5- spin electron , @xmath6 is an electron density operator , @xmath7 represents the hopping of electrons between site @xmath8 and @xmath9 , and @xmath10 are on - site ( off - site ) coulomb repulsion . 
   bracket @xmath11 denotes the next - neighbor pair , @xmath12 is the strength of the spin - orbit interaction and @xmath13 is the @xmath14=@xmath15-spin operator . in eq.([hso ] ) , the @xmath16-th site is in the middle of the next nearest neighboring pair @xmath8 and @xmath9 as shown in fig . [ 
   fig : honeycomb ] , and @xmath17 is the vector from the site @xmath8 to @xmath9 . 
   we start with the hopping matrix @xmath7 in eq.([eq : h0 ] ) for the bond connecting a pair of the nearest - neighbor sites @xmath18 , @xmath19 as the simplest extended hubbard model on the honeycomb lattice . 
   later , we will examine the effect of third neighbor hoppings @xmath20 . 
   we take @xmath21 as the unit of energy and set @xmath22 throughout the rest of this paper . for the non - interacting limit , @xmath23 , 
   the system becomes a topological insulator when @xmath12 is nonzero , which is identical to the topological phase of the kane - mele model  @xcite . 
   for the off - site coulomb repulsion , we mainly consider the second neighbor interaction ( @xmath24 ) , which is necessary for the emergence of the correlation - induced topological insulator  @xcite . 
   the second neighbor coulomb repulsions @xmath25 effectively generate the spin - orbit interactions , which are identical to @xmath12 , and induce topological insulator phases even for @xmath26 . 
   we note that the coulomb repulsion of on - site or the nearest neighbor site do not affect the stability of the tmi phase at the level of the mean - field approximation . 
   indeed , our mvmc results show that this is essentially true beyond the level of mean - field approximation which will be discussed in the later section . therefore , for the moment , we focus only on the effect of @xmath27 for the consideration of interaction effects .    . 
   next nearest neighbor interaction @xmath25 and the loop current @xmath28 flowing between the next nearest neighbor are shown by red and blue arrows . 
   , width=245 ]    here the tmi is the broken symmetry phase characterized by the order parameter @xmath28 defined by @xmath29 where the self - consistent mean fields for the second neighbor bonds are given by @xmath30 here , @xmath31 denotes the expectation value for next nearest neighbor bonds and the order parameter @xmath28 is physically interpreted as spin dependent loop currents flowing within a hexagons constituting the honeycomb lattice .    at the level of the mean - field approximation , 
   this quantum phase transition is understood by decoupling two - body electron correlation term of the next nearest neighbor bond , @xmath32 , as @xmath33 we also note that this phase transition to tmi , which is proposed not only on the honeycomb lattice but also in several other lattice models , belongs to an unconventional universality class , which depends on the dimension of the system and the dispersion of the electron band@xcite .      in this section 
   , we pursue the topological mott phase transition by employing the mean - field analysis and the variational monte carlo method . for the latter method 
   , we use a trial wave function of the gutzwiller - jastrow form , @xmath34 with a one body part , @xmath35 where @xmath36 is the variational parameters and @xmath37 is the number of the electrons in the system . 
   though this form of the wave function restricts itself to the hilbert subspace with the zero total @xmath38-component of @xmath14=@xmath15 , @xmath39 , it can describe topological phases on the honeycomb lattice as long as we use complex variables for @xmath36 . here , 
   @xmath40 and @xmath41 are the gutzwiller and jastrow factors defined as @xmath42 and @xmath43 respectively , with which the effects of electron correlations are taken into account beyond the level of mean - field approximation . 
   the expectation value , @xmath44 , is minimized with respect to variational parameters , @xmath36 , @xmath45 , and @xmath46 , by using the monte - carlo sampling and using the stochastic reconfiguration method by calculating gradient of the energy and the overlap matrix in the parameter space@xcite . 
   we optimize the parameters by typically 2000 stochastic reconfiguration steps .      in the present implementation of the variational monte carlo method with complex variables , 
   the feasible system size for the calculation is about up to 300 , from which we speculate properties in the thermodynamic limit . for this purpose 
   , we perform the size extrapolation using the standard formula @xmath47 where @xmath48 is the linear dimension of the system size . 
   we note that the order of taking the limit in the right hand side of eq.([zeta_limit ] ) is also important for the validity of the extrapolation , which is known as the textbook prescription for the defining spontaneous symmetry breakings . 
   another practical way to determine the spontaneous symmetry breakings is in principle the finite size scaling of the correlation for the order parameter . 
   however , the latter finite size scaling is not practically easy . 
   the correlation of @xmath28 becomes too small because @xmath28 itself is about the order of @xmath49 and the correlation becomes the order of @xmath50 which becomes comparable to the statistical error of the monte carlo sampling .    for the size extrapolation , we fit the data of the finite size calculations by a polynomial of @xmath51 , that is , we assume the size dependence by @xmath52 the above assumption for the finite size scaling is based on a practical observation and an analogy to the finite size scaling in the spin wave theory  @xcite . as a practical observation , @xmath53 for the limit @xmath54 and @xmath55 
   is scaled by @xmath51 . 
   we examine the ground states of the extended hubbard model on the honeycomb lattice by tuning the on - site coulomb repulsion @xmath56 and the second neighbor coulomb repulsion @xmath25 . even for the parameter sets favorable to the tmis that have been studied in the pioneering works on the tmi  @xcite 
   , we show that the tmi is not stabilized when we take into account other competing orders overlooked in the literature . 
   however , here , we reveal that , by tuning the fermi velocities of the dirac cones , the tmi is indeed stabilized .      in this section 
   , we consider long - period cdw states and show that 6-sublattice order is stable when @xmath25 becomes dominat . 
   this state is schematically shown in fig.[fig : cdw6sub](a ) , where the electron density per site is disproportionated into four inequivalent values . 
   when we pick up very rich sites or very poor sites , they constitute triangular lattices . 
   the moderately rich or poor sites constitute the honeycomb lattices . 
   sites labelled by 1 ( red ) , 2 ( light red ) , 3 ( light blue ) and 4 ( blue ) correspond to the sites where the electron densities are very rich , rich , poor and very poor , respectively . 
   ( b ) the growth of order parameter of the cdw state . 
   ( c ) resulting phase diagram of honeycomb lattice for @xmath56 and @xmath25 at the level of mean - field approximation . 
   we do not find the region of stable tmi phase . 
   , width=302 ]    figure [ fig : cdw6sub](b ) shows the growth of the order parameters of the 6-sublattice cdw state for several different parameters of @xmath56 . in this calculation , we have used the mean - field approximation , where we defined the mean field by @xmath57 where @xmath58 is defined through @xmath59 here , the condition of half filling is satisfied by the following identity , @xmath60    as can be seen from fig.[fig : cdw6sub](b ) , the order parameter grows above @xmath61 even if we set @xmath62 . 
   this is a quite serious problem for the stabilization of the tmi since the critical value of @xmath25 for the tmi at the mean - field approximation is also about @xmath63 . 
   in addition , the energy gain due to the formation of the topological mott order is much smaller than that for the cdw state . 
   when we consider larger values of @xmath56 , then the system becomes an antiferromagnetic state . 
   furthermore , the mean - field approximation often overestimates the ordered phase , because it does not take into account fluctuation effects seriously . actually , by using the mvmc method , we could not find the region where the topological insulator becomes the ground state in the parameter space of @xmath56 and @xmath27 . the resulting phase diagram by the mean - field approximation is shown in fig . 
   [ fig : cdw6sub ] ( c ) . there 
   the antiferromagnetic phase is denoted as nel , where the spin on different sites of the bipartite aligns in the anti - parallel direction . 
   as examined above , we found that the cdw state dominates and could not find parameter regions where the tmi is stable . 
   however , we find that the tmi becomes stable by modulating the fermi velocity at the dirac cones in the electronic band dispersion actually it is confirmed that , when the fermi velocity is 0 and then the dirac cones change to the quadratic band crossing points , the phase transition from a zero - gap semiconductors to a tmi occurs with infinitesimal coulomb repulsions  @xcite . 
   for the honeycomb lattice , it is possible to change dirac semimetals to quadratic band crossings by introducing the third neighbor hopping @xmath20 ( schematically shown in fig.[fig : t3 ] ( a ) ) . 
   then the part of the hamiltonian @xmath64 is replaced with @xmath65 where the third neighbor hoppings are given as @xmath66 here the tnn stands for the third neighbor . 
   we find that the fermi velocity linearly decreases by introducing @xmath20 and becomes @xmath67 at @xmath68 as @xmath69 , which is also shown in fig.[fig : t3 ] ( b ) . 
   though tuning nominal value of @xmath70 in the graphene is difficult , tuning of the fermi velocity at the dirac cone has been proposed in bilayer graphene by changing the relative orientation angle between two layers@xcite , which is effectively equivalent to the tuning of @xmath70 .     and 
   fermi velocity @xmath71 at dirac point . here 
   , @xmath71 is shown by the ratio to the fermi velocity at @xmath72 . , width=321 ]    figure [ fig : mfphase ] shows the phase diagram calculated by the mean - field approximation for @xmath73 , @xmath74 , and @xmath75 . 
   compared to the nel ordered phase and the tmi , the cdw phase is not largely affected by the fermi velocity . 
   therefore , the region of tmi recovers . 
   ( a ) 0.3 , ( b ) 0.4 and ( c ) 0.45 . 
   sm denotes the semimetal . 
   , width=321 ]      next , we examine the stability by using the mvmc method . in our mvmc method , we impose the translational symmetry on the variational parameters to reduce the computational cost , as @xmath76 here , @xmath77 is taken as any bravais vector of the 6 sublattice unit cell illustrated in fig.[fig : cdw6sub](a ) , in order to examine possible spontaneous symmetry breakings including the cdw shown in fig.[fig : cdw6sub](a ) and the tmi on an equal footing . 
   therefore the number of the sites for each calculation is taken as @xmath78 . 
   we have calculated for @xmath79 , where about 2500 variational parameters are used for the calculation of the largest size . to perform the extrapolation to the small external filed limit , @xmath80 , we have also calculated for several different strengths of the spin - orbit interaction , @xmath81 , and @xmath82 . 
   figure [ fig : vmcdata](a ) shows the numerical results for the order parameter of the tmi at @xmath83 and @xmath84 . 
   the sudden drops in @xmath28 around @xmath85 signals the emergence of the 6-sublattice cdw state . indeed , 
   when @xmath27 further increases , @xmath28 vanishes . 
   this is physically quite natural because @xmath28 is interpreted as the loop current and hard to be stabilized inside the cdw phase where electrons are locked at specific sites .     for @xmath86 at @xmath84 . 
   ( b ) size extrapolation at @xmath87 for several different values of @xmath27 . here 
   , four different sizes ( @xmath79 ) are calculated . 
   lines are results of quadratic fittings . 
   ( c)extrapolation of @xmath12 is shown for different values of @xmath27 . 
   lines are results of quadratic fitting . , 
   width=321 ]    same calculations are carried out for @xmath88 , which is employed in the size extrapolation of @xmath28 . 
   figure [ fig : vmcdata](b ) shows the size extrapolation at @xmath87 as a typical example . 
   since the cdw becomes dominant at @xmath85 , the results for @xmath89 are shown in fig.[fig : vmcdata](b ) . 
   when @xmath25 is small , @xmath28 empirically follows the size dependence @xmath90 with a constant @xmath91 as we see in fig.5(b ) . 
   this extrapolation is performed for four different values of @xmath12 , and then we get the values for @xmath92 . 
   these data are shown in fig . 
   [ fig : vmcdata](c ) , where the second extrapolation to @xmath54 is performed . 
   final results for @xmath93 are shown in fig . 
   [ fig : vmcresult](a ) , where the results of @xmath94 are also shown . 
   the errorbars are defined from the errors of the extrapolation , where the largest error of the first extrapolation is added to the errors of the second extrapolation . 
   we note that the error bars arising from the statistical errors of the monte carlo sampling are much smaller . 
   relatively large errorbars for small @xmath27 at @xmath94 is possiblly because of the existence of the critical point at a finite value of @xmath27 , which is about @xmath95 . 
   when @xmath27 exceeds this critical value and @xmath28 in the thermodynamic limit remains nonzero , the error becomes smaller as can be seen from fig . 
   [ fig : vmcresult](a ) . for @xmath83 , the order grows from small value of @xmath96 . 
   however , at @xmath94 , non - zero @xmath28 can not be detected for small values of @xmath27 . 
   though the estimate of the universality class from these data is difficult , theoretically it is expected to belong to that of the gross - neveu model@xcite , and our result does not contradict this criticality . for @xmath97 , we do not find the value of @xmath27 where @xmath28 in the thermodynamic limit remains nonzero , and therefore phase transitions is not expected . 
   this is shown in fig . 
   [ fig : vmcresult](b ) , where the size extrapolation at @xmath87 is shown . 
   there , @xmath28 becomes @xmath67 at @xmath98 for all @xmath27 , which is completely different from the behavior at @xmath94 and @xmath83 . the resulting phase diagram is shown in fig . 
   [ fig : ueffect](a ) . 
   although we expect that the mvmc results show larger critical values @xmath99 for the transition at @xmath100 , the estimated results indicate @xmath101 slightly smaller than the mean - field results shown in fig.[fig : mfphase ] . 
   the reason that @xmath101 becomes smaller at @xmath100 in the mvmc results is probably an artifact arising from a peculiar size dependence near the essential singularity at @xmath102 , as we see even in the mean - field calculation shown in fig . 
   [ fig : ueffect](b ) where , the size dependence in the mean - field calculation are shown . 
   the possible errors in in the estimate of @xmath28 is as large as @xmath103 and the resultant errors in the estimate of @xmath101 is around 0.1 . 
   therefore , the stability of the tmi phase over the cdw and nel phases in the region @xmath104 for @xmath100 is robust . 
   here we note that the boundary between the cdw and tmi phases around @xmath105 does not change when we take into account the quantum fluctuations ( by calculating with the vmc method ) , because it is a strong first - order transition . 
   are shown . 
   ( b)size extrapolation at @xmath106 and @xmath87 for several different values of @xmath27 . here , four different sizes ( @xmath79 ) are calculated . 
   lines are results of quadratic fittings . , width=321 ]    now we discuss the effect of the on - site coulomb repulsion . though the onsite coulomb repulsion does not affect the value of @xmath28 and therefore stability of the tmi in the mean - field approximation , our mvmc result shows that increasing @xmath56 decreases the value of @xmath28 if @xmath27 is fixed as shown in fig . 
   [ fig : ueffect](c ) . 
   while the increasing @xmath56 quantitatively decreases the value of @xmath28 , its effect does not destroy the stability of the tmi phase at @xmath107 at least if @xmath108 . 
   it also suppresses the emergence of the cdw . 
   therefore it may help to enlarge the region of the tmi phase . 
   the same effect is expected when we consider the coulomb repulsion for the nearest neighbor sites @xmath109 . 
   that is , it decreases the value of @xmath28 beyond the level of the mean - field approximation but does not essentially affect the phase transition . 
   however , we also note @xmath109 may cause another type of cdw , and stabilization of tmi should be examined against this cdw phase when @xmath109 is large .     and @xmath70 at @xmath110 obtained by mvmc calculations . 
   ( b)size dependence of mean - field calculation for several different values of @xmath27 at @xmath83 . 
   ( c ) @xmath27 dependence of order parameter of tmi for @xmath111 and @xmath112 at @xmath113 . 
   decrease in @xmath28 is observed by introducing @xmath56 . 
   , width=321 ]    as a qualitative difference from the mean - field approximation , we found that effect of @xmath56 decreases the value of the order parameter of the tmi . in the case of the border between the semimetal and the nel ordered phases , 
   it is expected that @xmath114 , namely the critical values for @xmath56 becomes larger by treating fluctuation effects carefully . 
   furthermore , we expect that @xmath114 becomes a function of @xmath27 , which may enhance the fluctuation of the nel order and suppress the phase transition . in our calculation 
   , we did not find nel ordered states or mixed states of tmi and nel ordered states . on the other hand , 
   the cdw phase is expected to be much more stable against fluctuations and the mean - field solution gives a reasonably good description because of a large scale of the energy gain for the cdw phase in comparison to other phases as adequately shown even by the mean - field approximation . at the boundary of the tmi to nel 
   ordered phases , the universality class may change as suggested in the kane - mele - hubbard model@xcite . 
   here , we discuss the realization of the tmi in the real solids . a primary candidate of tmis is graphene . as a well - known fact , 
   graphene is nothing but a two - dimensional honeycomb network of carbon atoms and hosts dirac electrons . 
   however , it is also a well - known fact that , in free - standing graphene and graphene on substrates , significant single - particle excitation gaps have not been observed yet@xcite . 
   below , we examine possible routes toward the realization of the tmi in graphene - related systems . 
   first of all , as already studied above , the suppression of the fermi velocity of the dirac electrons is crucial for the stabilization of the tmis . 
   as extensively studied in the literature  @xcite , twisted bilayer graphene ( tblg ) offers dirac electrons with tunable fermi velocities . by choosing stacking procedures , the quadratic band crossing , in other words , the zero fermi velocity limit , 
   is also achieved , which has been already observed in experiments@xcite . 
   next , we need to clarify competitions with any other possible symmetry breakings in the tblg with the fermi velocity smaller than that of graphene . 
   the suppressed fermi velocities may possibly cause instabilities towards not only the tmis but also other competing orders as discussed in this paper . 
   for clarification of the competition we need an _ ab initio _ estimation of effective coulomb repulsions which directly correspond to the coulomb repulsions in the extended hubbard model@xcite . 
   the _ ab initio 
   _ study on the effective coulomb repulsions employs a many - body perturbation scheme called constrained random phase approximation ( crpa)@xcite . the crpa estimation gives the following values for the coulomb repulsions : the on - site and off - site coulomb repulsions are given as @xmath115 , @xmath116 , and @xmath117 with @xmath118 ev for free - standing graphene . 
   if we neglect longer - ranged coulomb repulsions , we expect the nel or cdw orders by employing these crpa estimates of @xmath119 , and @xmath25 in graphene and tblg . 
   therefore , the free - standing graphene and tblg do not offer a suitable platform for the tmis . 
   however , by choosing dielectric substrates , the strength of the coulomb repulsions , @xmath120 , and @xmath121 , is suppressed due to enhancement of dielectric constant as @xmath122 , and @xmath123 , where @xmath124 is defined by dielectric constants of each materials as @xmath125 . 
   ( here , we ignored the possible reduction of the effective dielectric constants at small distances . ) then , we may approach the parameter region , where tmis become stable , as shown in fig.[fig : ueffect](b ) . 
   if we neglect @xmath126 , dielectric substrates with @xmath127 are enough to stabilize the tmis . 
   even when the nearest - neighbor coulomb repulsion @xmath126 is taken into account , tmi is expected to be stable as long as @xmath126 is not strong .    in the above discussion , we neglected further neighbor coulomb repulsions , namely , the third neighbor coulomb repulsions @xmath128 , the fourth neighbor ones @xmath129 , and so on . to justify the above discussion , we need to screen the further neighbor coulomb repulsions by utilizing dielectric responses of the substrates and/or adatoms . here 
   we note that the screening from atomic orbitals on the same and neighboring sites effectively reduces the coulomb repulsions as extensively studied by using crpa@xcite . 
   the relative strength of the on - site and second neighbor coulomb repulsions , @xmath130 , may also be controllable by utilizing the adatoms , which are expected to efficiently screen the on - site coulomb repulsions if the adatoms is located just on top of the carbon atoms . 
   if we combine the control of the @xmath130 with the suppression of the further - neighbor coulomb repulsions , the above discussion may be relevant . by utilizing the adatoms , 
   the nearest - neighbor coulomb repulsions @xmath126 are also expected to be well - screened by adatoms on nearest - neighbor carbon - carbon bonds . 
   the suppression of @xmath126 is helpfull for suppressing the cdws competing with the tmi .    finally , we estimate the single - particle excitations gap @xmath131 induced by the tmi . 
   the excitation gap is crucial for actual applications of the tmi as a spintronics platform . 
   if we set @xmath132 and expect the tblg with @xmath133 ev , we obtain the excitation gap up to 0.1 ev , where we use the mean - field estimation of the gap @xmath134@xcite with the value of the mvmc result for @xmath28 and assumed this formula is valid in the presence of electron correlation . 
   the estimated gap scale is substantially larger than the room temperatures . 
   our theoretical results support that topological insulators with such a large excitation gap @xmath135 ev are possibly obtained by using abundant carbon atoms . 
   in this paper , we have studied the realization of the tmi phase for the electronic systems on a honeycomb lattice by using the mean - field calculation and mvmc method . 
   we found that the cdw of the 6 sublattice unit cell is much more stable than the previously estimated cdw with smaller unit cells for the simplest case where the electronic transfer is limited to the nearest neighbor pair . 
   for the stabilization of the tmi we need to suppress the fermi velocity at the dirac point than the standard dirac dispersion for the case with only the nearest neighbor transfer . in the case of the honeycomb lattice , this is realized by introducing the third neighbor hopping @xmath20 and we have given quantitative criteria for the emergence of the tmi . 
   related real material is a bilayer graphene where the fermi velocity is tuned by changing the rotation angle between two parallel layers@xcite . 
   actually , the quadratic band crossing is realized when the rotation angle is @xmath67 ( known as the ab stacking bilayer graphene@xcite ) , which is mimicked by @xmath136 . since smaller values of fermi velocity stabilizes the tmi at smaller values of @xmath27 , 
   its effective control may offer a breakthrough in the realization of two dimensional tmis . 
   we need further analyses for experimental methods of controlling the stability of the tmis and _ ab initio _ quantitative estimates of the stability for bilayer graphenes , which are intriguing future subjects of our study . 
   the authors thank financial support by grant - in - aid for scientific research ( no . 
   22340090 ) , from mext , japan . 
   the authors thank t. misawa and d. tahara for fruitful discussions . a part of this research 
   was supported by the strategic programs for innovative research ( spire ) , mext ( grant number hp130007 and hp140215 ) , and the computational materials science initiative ( cmsi ) , japan .
summary: realization and design of topological insulators emerging from electron correlations , called topological mott insulators ( tmis ) , is pursued by using mean - field approximations as well as multi - variable variational monte carlo ( mvmc ) methods for dirac electrons on honeycomb lattices . 
   the topological insulator phases predicted in the previous studies by the mean - field approximation for an extended hubbard model on the honeycomb lattice turn out to disappear , when we consider the possibility of a long - period charge - density - wave ( cdw ) order taking over the tmi phase . 
   nevertheless , we further show that the tmi phase is still stabilized when we are able to tune the fermi velocity of the dirac point of the electron band . beyond the limitation of the mean - field calculation 
   , we apply the newly developed mvmc to make accurate predictions after including the many - body and quantum fluctuations . by taking the extrapolation to the thermodynamic and weak external field limit 
   , we present realistic criteria for the emergence of the topological insulator caused by the electron correlations . by suppressing the fermi velocity to a tenth of that of the original honeycomb lattice , 
   the topological insulator emerges in an extended region as a spontaneous symmetry breaking surviving competitions with other orders . 
   we discuss experimental ways to realize it in a bilayer graphenesystem .