SciPost Submission Page
Chiral pwave superconductivity in twisted bilayer graphene from dynamical mean field theory
by B. Pahlevanzadeh, P. Sahebsara, D. Sénéchal
This is not the latest submitted version.
This Submission thread is now published as
Submission summary
Authors (as registered SciPost users):  Peyman Sahebsara · David Sénéchal 
Submission information  

Preprint Link:  scipost_202103_00017v1 (pdf) 
Date submitted:  20210316 17:47 
Submitted by:  Sénéchal, David 
Submitted to:  SciPost Physics 
Ontological classification  

Academic field:  Physics 
Specialties: 

Approaches:  Theoretical, Computational 
Abstract
We apply cluster dynamical mean field theory with an exactdiagonalization impurity solver to a Hubbard model for magicangle twisted bilayer graphene, built on the tightbinding model proposed by Kang and Vafek [2018]. This model applies to the magic angle $1.30^\circ$. We find that triplet superconductivity with $p+ip$ symmetry is stabilized by CDMFT, contrary to other triplet or singlet order parameters. A minimum of the order parameter exists close to quarterfilling and threequarter filling, as observed in experiments.
Current status:
Reports on this Submission
Anonymous Report 3 on 2021429 (Invited Report)
 Cite as: Anonymous, Report on arXiv:scipost_202103_00017v1, delivered 20210429, doi: 10.21468/SciPost.Report.2855
Strengths
1 Very active subject in condensed matter Physics of 2D materials, with many open questions.
2 Complex numerical method that must be able to give meaningful results to properly study the electronic correlations.
3 The results shown seem to be in fairly good agreement with experimental results.
4 The manuscript is well written and seems understandable to me.
Weaknesses
1 The model studied is necessarily approximate, and it is difficult to know whether the results obtained depend much or not on this model.
2 Weak comparison of the results with recent literature.
Report
This manuscript proposes a numerical study of correlations and superconductivity in twisted bilayre graphene (TBG) with a rotation angle close to the first magic angle (MATBG). The results shown seem new and in good agreement with experimental results. A main result is the prediction of a stable superconducting states for ¼ filling and ¾ filling. The method used (cluster dynamical mean field theory (CDMFT)) should be effective in dealing with the electronelectron interactions that govern the physics of these systems. Thus, this work done with a complex numerical method brings new results, which confirm other simpler approaches. Therefore, I think this work deserves to be published provided that the authors answer the following requests.
Requested changes
1 CDMFT is difficult to implement numerically and it does not allow the study of large systems. This is why the authors have used a simplified version of the tightbinding Wannier model proposed by Ref. [1] to simulate the 4 lowenergy flat bands in MATBG. This way of doing is quite legitimate in view of the difficulty of the numerical method, but it should be shown that the simplifications made do not change too much the properties of flat bands. For instance, the authors neglect the hopping terms beyond the first neighbors. This approximation must be discussed further. A suggestion could be to show the effect of this approximation on noninteracting bands.
2 A number of recent results from the literature are cited in the introduction, but the results obtained are not really compared with them. It is not possible to compare with all that has been published on the subject, but it is important to situate the results in relation to works already published. The authors should improve this aspect of their manuscript.
3 A conclusion insisting on the originality and new results for understanding the electronic properties of flat bands in MATBG would help the reader.
Anonymous Report 2 on 2021418 (Invited Report)
 Cite as: Anonymous, Report on arXiv:scipost_202103_00017v1, delivered 20210418, doi: 10.21468/SciPost.Report.2803
Strengths
1) study of a timely subject motivated by experiment
2) the applied EDCDMFT with elaborate supercluster setup is stateoftheart
3) numerical instabilities and limitations of the technique are discussed
Weaknesses
1) discussion of the results with respect to experiments and numerical studies could be improved
2) some technical details should be added
Report
The manuscript presents a numerical analysis of an effective Hubbard model of twisted bilayer graphene (TBG) obtained from maximally localized Wannier functions. The method of choice, exact diagonalization based cluster dynamical meanfield theory (EDCDMFT), is well suited to study the question of superconducting instabilities of different symmetry at zero temperature as a function of electronic filling.
The authors have studied a minimal version of the Wannier model introduced in Ref.[1] by employing a supercluster setup consisting of four symmetryrelated copies of a 4site cluster that is embedded in an electronic bath represented by 6 discretized bath sites. Symmetrybreaking terms were added to the bath to test the model for existence of different singlet and triplet superconducting phases. After employing an elaborate procedure to ensure stability of the CDMFT selfconsistency cycle, the authors obtain a stable superconducting solution with p+ip symmetry away from halffilling, which is (in accordance with experiments on TBG) suppressed around 1/4 and 3/4 filling.
The paper is overall well written and includes some precious technical details on how to implement a stable CDMFT procedure that can describe triplet superconductivity for the presented model. In this respect, in particular the choice of clusterbath geometry based on symmetry considerations and the idea of using a counterterm in the impurity Hamiltonian to enforce spatial isotropy are noteworthy.
However, some technical aspects of EDCDMFT are not explained or referenced properly and some choices both on the model level and on the technique could be challenged. This concerns in particular the choice of hopping terms and the value of the interaction strength of their model in view of TBG as well as consequences of the specific clusterbath geometry and of neglecting spatially anisotropic solutions in the CDMFT procedure.
Before I can recommend the paper for publication in SciPost, I would therefore like to ask the authors to comment on the points listed below.
Requested changes
1) Model
1a  The authors consider a minimal tightbinding model, which they base on the maximally localized Wannier orbitals of Ref. [1]. Whereas this is helpful from a pedagogical point of view, restricting the tightbinding model to shortrange hopping terms can lead to some serious deviations from the DFT band structure of TBG, see e.g. Fig.5 of the similar Wannier study of Koshino et al [Koshino et al. PRX 8, 031087 (2018)]. If there exists no fundamental issue with including more terms of the Wannier model in their CDMFT study, the authors might want to comment on why they restricted the hopping terms and how they made the precise choice of hopping terms used in the manuscript. Most notably, Table 1 of the manuscript refers to them as being 'the most important amplitudes computed in Ref.[1]'. However, according to table I of Ref.[1] there are other hopping terms with comparable hopping strength and/or larger multiplicity (e.g. $t_{14}[1,0]=\dots=0.046+0.083i$).
1b  The choice of the local Hubbard interaction strength U requires some further explanation: In Fig.4 the authors use values of U ranging from 0.5 to 3.0 meV. This should be discussed in context of estimates from experiment (U larger than 30 meV possible for an angle of 1.3 degree in [Cao et al. Nature 556, 80 (2018)]) and from theory, e.g. obtained from a continuum model [Pizarro et al. PRB 100, 161102(R) (2019)] or from constrained random phase approximation [Goodwin et al. PRB 100, 235424 (2019); Vanahal and Pollet, PRB 102, 035154 (2020)].
2) Technique and results
2a  A few technical aspects of the technique should be mentioned or explained in more detail:
 The restricted Nambu formalism for the sc pairing fields in context of CDMFT should be either written down explicitly or referenced properly.
 Which is the typical number of iterations until convergence of the CDMFT selfconsistency cycle and how does the computation time of step (2) compare to that of step (4)?
 Which algorithm is used to optimize the set of bath parameters that minimizes the distance function (9) and does it preserve ergodicity? In other words, can the authors be sure that the final solution in the 13dimensional variational parameter space is not a local, but a global minimum of the distance function?
 For less complex clusterbath geometries, the convergence of EDCDMFT with respect to the number of bath sites has been already studied in detail. Here, only two bath sites per cluster site (at the boundary) are used. How can the authors be sure that this discretization is sufficient for a good representation of the hybridization function? Have they compared the quality of the CDMFT solution to the one obtained with three bath sites per cluster site? Such a geometry would amount to an effective 13site problem for the ED solver, which could still be feasible, even when treating symmetrybreaking terms like those of eqs (1) and (2).
2b  the authors might want to highlight the consequences of not including neither charge density modulations (CDW phases), nor antiferromagnetic solutions in the present study. In particular the latter could be rather easily included via a staggered Weiss field on the bath sites. Whereas the simultaneous treatment of SC and AFM would of course increase the already impressive variational parameter space and might thereby go beyond the scope of this paper, the study of the existence of a (pure) AFM phase should be feasible. The presence of AFM order at halffilling would provide a more complete phase diagram and underline the qualitative agreement of the studied model with experiments of TBG. Perhaps the authors could also comment on the possibility of further symmetrybreaking phases that might compete at n=1/4 and n=3/4 with the p+ip sc solution.
2c  the authors mention a weak singlet phase of (d+id) symmetry when not imposing spatial uniformity. However, they do not discuss this phase, arguing that it disappears when imposing the counterterm $H_{ct}$. This is very unfortunate, since one of the important discussions in the field is about this very question. CDMFT is one of the very few techniques that can not only treat sc phases of d+id or p+ip symmetry, but even allows to study the interplay of different symmetrybreaking instabilities. Since CDMFT finds both d+id and p+ip phases, a more thorough discussion of the d+id phase would be desirable.
3) Minor comments on spelling and formatting:
 typo on page 7: 'bybridization'
 perhaps the standalone Figure 1 could be integrated in the text
 semicolons (';') are a bit overused throughout the paper. This is in particular true for the last paragraph (conclusions), which needs proper rewriting.
Author: David Sénéchal on 20210612 [id 1502]
(in reply to Report 2 on 20210418)
We thank the referee for this very thorough review of our work.
1A. The referee suggests that we refine the tightbinding model by adding another set of hopping parameters that we neglected since we elected to keep the largest two sets of couplings on any given layer. We followed that suggestion and redid the computations including the next largest hopping amplitude. This explains in part why our response was delayed. The results are qualitatively the same for the triplet case but resurrect the singlet solution, which is now stabilized, although with a smaller order parameter and higher energy than the triplet solution. These additional hopping terms naturally lead to a minimum at fillings 0.5 and 1.5 without need for cluster counterterms, so that new technique was dropped from the paper.
1B. The referee asks us to justify the range of values of U used in the computations. We did not use any estimate from the references given by the referee to select our range of U (from 0.5 meV to 3 meV). We simply chose a range that would take us from intermediate coupling to strong coupling. In the revised version we situate values suggested by those references in that range. In particular, the value of U quoted by Cao et al (U~20 meV) seems very large in this context. We added the value U=5meV. This is near the lower bound of the estimates of Goodwin et al. Yet we are of the opinion that this value being already larger than the bandwidth of the model, no new physics will emerge by pushing U beyond this value.
2A. We bring precisions about the methods in the revised version, as suggested by the referee. However, we do not go as far as performing a new set of computations with 3 bath sites per boundary cluster site, which would bring the total number of orbitals to 13. This would take considerable resources given that (i) the cluster has no point group symmetry to be exploited easily (ii) complex numbers are needed, (iii) the exact diagonaliations have to be performed very many times. Also, our reasoning for the bath parameterization is inspired considering the bath sites as a sort of mockup of the neighboring sites, which leads to this impurity model of 6 bath sites per cluster.
2B. We agree with the referee that these additional orders are interesting to study and we indeed plan on this. In the context of the present paper, this would take us a bit too far. We are pretty much convinced that such a model would lead to antiferromagnetism at halffilling and we are currently investigating, in a different project, the possibility of charge order near quarter filling. But we use slightly different methods for this and including these considerations in the present paper would unduly delay publication as well as break the unity of the paper in terms of methods and impurity model.
2C. We agree. It turns out that including one additional set of subleading hopping terms has stabilized the singlet solution (d+id) and we cover this in the revised version.
Anonymous Report 1 on 2021412 (Invited Report)
 Cite as: Anonymous, Report on arXiv:scipost_202103_00017v1, delivered 20210412, doi: 10.21468/SciPost.Report.2780
Strengths
DMFT study of superconductivity in twisted bilayer Graphene
Weaknesses
 The effective model used in the manuscript cannot completely describe the lowenergy states in twisted bilayer Graphene.
 The results are poorly explained.
Report
The authors of the manuscript "Chiral pwave superconductivity in twisted bilayer graphene from dynamical mean field theory" study superconductivity in an effective model for flat bands emerging in twisted bilayer Graphene.
For this purpose, they use cluster DMFT combined with exact diagonalization, where the authors construct different impurity models depending on the order parameter studied. The authors find that chiral pwave superconductivity is the stable pairing in their calculations.
Currently, I cannot recommend this manuscript for publication in SciPost Physics for the following reason:
Most of the manuscript is about the method. The actual results section only shows a single figure, including the pwave order parameter for different filling and interaction strengths. There is not much interpretation or explanation of their results. For example, the authors say, "We expect the superconducting order parameter to be minimum, if not zero, at quarter filling (n=0.5), as observed in experiments.", but a theoretical explanation is not given.
Furthermore, superconductivity has been studied in effective models of twisted bilayer Graphene with different methods, and different results have been found, e.g., extended swave, pwave, and dwave superconductivity.
The authors should compare their results to these previous studies and discuss and explain differences and similarities. Why do the current results favor chiral pwave superconductivity?
Besides this, I have the following comments:
(1) Concerning the effective model: The authors use the effective 4band model proposed in Phys. Rev. X 8 031088. Recently, the importance of fragile topology for the flat bands in twisted bilayer Graphene has been discussed. For example, it has been discussed that there cannot exist an effective 4band tightbinding model for the flat bands fulfilling all symmetries (see, e.g., Phys. Rev. B 99, 195455 ) due to the topology in the bands. Thus, the current model breaks certain symmetries present in the real material. I believe this should be discussed shortly in the manuscript.
Furthermore, the effective basis (Wannier orbitals) in the current model reaches over many atoms in the real material. How realistic is the chosen local interaction and its strength?
(3) Concerning the method: I would appreciate it if the authors could give more information about the actual exact diagonalization. What is the dimension of the Fock space? What symmetries have been incorporated into the Hamiltonian? How was the spectral function calculated?
(3) As mentioned above, the authors should better explain their results.
Why has the order parameter a minimum at n=0.5?
Why does the order parameter become constant for small fillings?
For U=3meV, the order parameter vanishes for 0.55<n<0.7; why?
Is there an explanation why other pairings, such as dwave, vanish in the calculations?
Author: David Sénéchal on 20210612 [id 1501]
(in reply to Report 1 on 20210412)

The main issue raised here is the apparent contradiction between the use of a fourband tightbinding model proposed in Phys. Rev. X 8 031088 and a argument given in Phys. Rev. B 99, 195455 that a fourband model cannot sustain the fragile topology expected in the lowenergy model. Please see the discussion to that effect in Phys. Rev. B 102, 075142. While it is true that using localized Wannier functions with a purely local action of the symmetries is strictly incompatible in the continuum limit with the fragile topology, the error committed doing this is small. In any case, this issue is dwarfed by the approximations we make by keeping only the dominant hopping terms. In addition, we expect that introducing strong interactions would destroy any fragile topology that would be preserved in a noninteracting model. Hence, for the purpose of identifying the most favored superconducting pairing in a strongly interacting model, we believe this should not be an issue.

A concern is the use of local interactions when in reality we expect the interaction to extend over neighboring Wannier states. This a very legitimate concern and is the main weakness of our work. It is rooted in the small size of the cluster used in our CDMFT computations. We could keep extended interactions contained in a cluster, but this would ignore intercluster extended interactions. A way to alleviate this problem would be to implement a form of intercluster meanfield theory. We mention this possibility in the revised manuscript, but reserve its implementation for future work.

The referee asks for some precisions about the methodology: the exact diagonalization, the spectral function, etc. We provide some of this in the revised version.

The referee wonders why we expect the order parameter to have a minimum at the quarter filling (n=0.5). We expect this because (i) this is observed in experiments at magic angle 1.08 degrees and is the consequence of an insulating state at that commensurate filling value. Commensurate filling is more likely to lead to an insulating state because of interactions, in particular extended interactions. In our study, even if extended interactions are not taken into account a sharp minimum in the superconducting order parameter is observed. Note that our revised results no longer use counterterms to enforce the precise location of this minimum: the additional hopping terms included suffice.
Author: David Sénéchal on 20210612 [id 1503]
(in reply to Report 3 on 20210429)The model we use does not exactly neglect hopping beyond first neighbors. In the version of the manuscript reviewed by the referees, hopping was kept with first and third neighbors. In the revised version we have added hopping with fourth neighbors. Third and fourth neighbor hopping are not present within the small foursite cluster but still have an effect through the selfconsistency relation. To clarify this point, we have added a figure that sketches the hopping terms used in the paper. We also added a remark to the effect that hopping terms beyond nearest neighbor have an effect through the CDMFT selfconsistency relation.
The revised version of the manuscript comments on the relations of our results with other published results for the superconducting order parameter.
The conclusion has been expanded.