{'fraction_non_alphanumeric': 0.03582137936623713, 'fraction_numerical': 0.008104384471999352, 'mean_word_length': 4.458748064587481, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generalized Pomeranchuk instabilities in graphene

We study the presence of Pomeranchuk instabilities induced by interactions on a Fermi liquid description of a graphene layer. Using a recently developed generalization of Pomeranchuk method we present a phase diagram in the space of fillings versus on-site and nearest neighbors interactions. Interestingly, we find that for both interactions being repulsive an instability region exists near the Van Hove filling, in agreement with earlier theoretical work. In contrast, near half filling, the Fermi liquid behavior appears to be stable, in agreement with theoretical results and experimental findings using ARPES. The method allows for a description of the complete phase diagram for arbitrary filling.


INTRODUCTION
Correlated electron systems in two dimensions have attracted a lot of attention in the last years, especially due to an important number of experiments that provide undisputable evidence of the existence of new exotic phases of matter.
One such example corresponds to nematic and stripe (smectic) phases in high T c superconductors in the underdoped region and fractional quantum Hall effect systems at high magnetic fields [1]. A nematic phase is characterized by orientational but not positional order and it has been proposed to explain the observed transport anisotropies. One important point about these phases is that they arise spontaneously, decreasing the rotational symmetry without a lowering of the lattice symmetry. Another more recent case is given by strontium ruthenate Sr 3 Ru 2 O 7 , which is well modeled as a bilayer system and shows a large magnetoresistive anisotropy. This observation has been argued to be consistent with an electronic nematic fluid phase. Experimentally, two consecutive metamagnetic transitions have been observed and the region in between has been proposed to be a consequence of a Pomeranchuk instability, due to a nematic deformation of the Fermi surface, in very close analogy to what happens in fractional quantum Hall gallium arsenide systems [2]. Yet another interesting material is the heavy fermion compound URu 2 Si 2 in which a hidden order phase arises through a second order transition at around 17.5 0 K. The order parameter of this new phase has remained elusive to theorists up to date. Different types of order have been proposed, but the situation is still under debate [3]. In recent work, Varma and Zhu [4] have proposed that this phase transition could correspond to a Pomeranchuk instability inducing a deformation in the antisymmetric spin channel, stabilized by a phase characterized by a helicity order parameter.
The experimental findings mentioned above triggered different theoretical studies on low dimensional correlated systems to search for such exotic phases [1]. Special attention has been paid to the possibility of Pomeranchuk instabilities [5] giving rise to such novel phases [6,7,8,9,10,11,12,13,14,15,16]. In a previous paper [17], motivated by these investigations, we developed a generalization of Pomeranchuk's method to search for instabilities of a Fermi liquid. The method we presented is applicable to any two dimensional lattice model with an arbitrary shape of the Fermi surface (FS) at zero temperature. The main results of our previous paper were summarized in the form of a recipe whose steps we give below for completeness. Our method is particularly well suited to analyze systems with weak interactions and then graphene appears as a perfect arena to test our techniques, since electron interactions are argued to be small due to strong screening.
The recent isolation of graphene [18], the first purely two-dimensional material, which is made out of carbon atoms arranged in a hexagonal structure, led to an enormous interest and a large amount of activity in studying its properties. The low doping region, near to half filling, became the subject of attention due to the peculiar behavior described by chiral massless charge carriers. Several Van Hove singularities are present at energies of the order of the hopping parameter E ∼ 2.7 eV and these singularities are expected to have an important role in the properties of the system. Although in first approximation graphene layers are well modeled by free fermions hopping on a hexagonal lattice, there have been a number of works in the literature where the effects of electron-electron interactions were taken into account [19,20]. However, such analysis were centered in the undoped (half-filling) regime or very close to it, and explicit analytic results for a wide range of fillings are still lacking.
In the present paper we apply the method developed in [17] to fermions in a graphene sheet in the presence of electron-electron interactions up to nearest neighbors. A previous work using a mean field approach, showed that long range interactions could lead to Pomeranchuk instabilities [21] at the Van Hove filling. On the other hand, at dopings very close to half-filling, it was anticipated that graphene should behave as a Fermi liquid [22]. This was later confirmed experimentally by ARPES exploration of the FS [23]. The results presented here are consistent with these findings, while our method allows for a more complete and systematic study of the whole space of fillings.

TWO DIMENSIONAL POMERANCHUK INSTABILITIES: REVIEW AND IMPROVEMENT OF THE METHOD
In this section we will review the generalization of Pomeranchuk method first presented in [17], and discuss a shortcut that can be used as an alternative to the change of variables proposed there.

Review of the method
According to Landau's theory of the Fermi liquid, the free-energy E as a functional of the change δn k in the equilibrium distribution function at finite chemical potential µ can be written, to first order in the interaction, as Here ε(k) is the dispersion relation that controls the free dynamics of the system, the interaction function f (k, k ′ ) can be related to the low energy limit of the two particle vertex. Note that we are omitting spin indices, and considering only variations of the total number of particles n k = n ↑k + n ↓k . This implies that the considerations that follow will be valid in the absence of any external magnetic field and at constant total magnetization. The identification of these two functions is the starting step of our calculation: Step 1: write the energy as in (1) and identify the functions ε(k) and f (k, k ′ ).
Instead of the cartesian variables in momentum space (k x , k y ) we will find convenient to define a new set of curvilinear variables (g, s) according to The variable g varies in the direction normal to the Fermi Surface (FS), whose position is defined by g = 0. We choose s such that it is constant at constant distance to the FS, varying in the longitudinal direction tangent to the FS, namely it satisfies the restriction ∇s(k).∇g(k) = 0 .
Since solving this for s may be a difficult task, we develop bellow an alternative procedure that replaces this calculation. Even if not mandatory, we will chose the variable s such that we give a complete turn around each connected piece of the FS when it runs from −π to π. From the above change of variables we obtain the Jacobian which is the relevant outcome of this step of the calculation: Step 2: with the help of the dispersion relation ε(k) obtained in step 1, change the variables according to (2) to obtain the Jacobian (4).
In a stable system, the energy (1) should be positive for all δn k that satisfy the constraint imposed by Luttinger theorem [24] d 2 k δn k = 0 .
Pomeranchuk's method roughly consists on exploring the space of solutions of constraint (5) to find a δn k that turns the energy into negative values, thus pointing to an instability of the system.
In terms of our new variables g and s, we can write δn k(g,s) at zero temperature as where H is the unit step function and δg(g, s) is an small perturbation parameterizing the deformation of the FS. Replacing into the constraint (5), changing the variables of integration according to (2) and performing the integral to lowest order in δg we get dsJ(s)δg(s) = 0 .
Here J(s) = J(g, s)| g=0 and δg(s) = δg(g, s)| g=0 . In case the FS has a nontrivial topology, the integral includes a sum over all different connected pieces. We see that in order to solve the constraint δg(s) can be written as in terms of a free slowly varying function λ(s). Even if in (5) and (7) a sum over different connected pieces of the FS may be assumed, this particular solution does not consider excitations in which some particles jump between different connected pieces .
Using the change of variables (2) and with the help of eqs. (6) and (8), we can rewrite the energy E to lowest order in δg as where we call ∂ s λ(s) = ψ(s) and f (s, Note that the left hand side of the stability condition E > 0 has two terms, the first of which contains the information about the form of the FS via J −1 (s), while the second encodes the specific form of the interaction in f (s, s ′ ). There is a clear competition between the interaction function in the second term of the integrand and the first term that only depends on the geometry of the unperturbed FS.
We see in (9) that the energy E is a bilinear form acting on the real functions ψ(s) that parameterize the deformations of the FS This is the next step of our calculation that can be summarized as: Step 3: write the energy as the bilinear form (9) using the functions J(s) and f (s, s) identified in steps 1 and 2.
The stability condition is then equivalent to demanding this bilinear form to be positive definite for any possible smooth deformation of the FS where L 2 [FS] stands for the space of square-integrable functions defined on the Fermi surface. In consequence, a straightforward way to diagnose an instability is to diagonalize this bilinear form looking for negative eigenvalues. To see that, we choose an arbitrary basis of functions {ψ n } n∈N of L 2 [FS], in terms of which we can write ψ(s) = n a n ψ n (s) , (12) and the the stability condition now reads This configures our Step 4: choose an arbitrary basis {ψ n } n∈N of the space of functions L 2 [FS].
The above defined bilinear form can be considered as a pseudo-scalar product in L 2 [FS]. In general the functions of the basis {ψ n } n∈N will not be mutually orthogonal with respect to this product. Moreover in the presence of instabilities, the pseudo-scalar product may lead to negative pseudo-norms ψ, ψ < 0.
We can then make use of the Gram-Schmidt orthogonalization procedure to obtain a new basis of mutually orthogonal functions {ξ n } n∈N . In terms of them an arbitrary deformation of the FS parameterized by a function ψ(s) can be decomposed as which implies that the stability condition on such deformation will read In summary: Step 5: apply the Gram-Schmidt orthogonalization procedure to go from the arbitrary basis {ψ n } n∈N chosen on step 3 onto a basis of mutually orthogonal functions {ξ n } n∈N .
In (15) we note that the only possible source of a negative sign is in the pseudo-norms χ n = ξ n , ξ n . In case the i-th pseudo-norm χ i is negative, a deformation parameterized by ψ(s) ∝ ξ i (s) is unstable. In other words the pseudo-norms {χ n } n∈N can be taken as the stability parameters, a negative value of χ i indicating a instability in the i-th channel. Then: Step 6: compute the pseudo-norms of the new basis functions {ξ n } n∈N . If for a given channel i one finds that These six steps summarize the generalized Pomeranchuk method. It can be applied to any two dimensional model with arbitrary dispersion relation and interaction. Note that since L 2 [FS] is infinite dimensional, the present method is not efficient to verify stability: at any step i it may always be the case that for some j, χ i+j < 0. Moreover, in the case of nontrivial topology, excitations consisting on particles jumping between different connected pieces of the FS were not included in the solution of the constraint (7), and they may lead to additional instabilities. The same is true excitations involving spin or color flips, that we are not considering.

An alternative to the change of variables
As advanced, to avoid the task of solving the constraint (3) that defines the variable s, we will develop here an alternative procedure to derive the form of the Jacobian evaluated on the FS J(s).
We begin by defining a parametrization of the FS in terms of an arbitrary parameter t. In other words, given the function g(k) defined in (2), we choose k(t) such that ∀t : g(k(t)) = 0. In terms of such parametrization we can decompose the Dirac δ function as (a proof of this formula is given in the Appendix).
The integral I of an arbitrary function F (g, s) along the FS can be written as or in other words Changing variables to k I = d 2 k F (s(k), g(k)) J −1 (s(k), g(k)) δ(g(k)) , (20) and replacing (17) we get or, interchanging the order of the integration and performing the k integrals I = dt F (s(k(t)), g(k(t))) J −1 (s(k(t)), g(k(t))) |k(t)| |∇g(k(t))| .
Now using the fact the g(k(t)) = 0 and defining the parameter t such that s(k(t)) = t, we get that in order to be compatible with (18) imply Then with all the above, we can replace our previous step 2 by a new version Step 2': with the help of the dispersion relation ε(k) obtained in step 1, parameterize the FS and obtain the Jacobian evaluated on the FS according to (24).
Then, even if it may be very difficult to solve the partial differential equation (3) in order to explicitly obtain the Jacobian, its restriction to the FS is all what we need, and can be obtained by the much easier task of parameterizing the FS.

POMERANCHUK INSTABILITY IN GRAPHENE
In the present section we will apply the method reviewed in the previous pages to the case of fermions in a graphene layer with Coulomb interactions.

Free Hamiltonian: tight-binding approach
Graphene is made out of carbon atoms arranged in a hexagonal lattice. It is not a Bravais lattice but can be seen as a triangular lattice with two atoms per unit cell. The tight-binding Hamiltonian for electrons in graphene considering that electrons can hop only to nearest neighbor atoms has the form where a σi , b σi are the creation and annihilation operators related to each of the unit cell atoms. The diagonalized Hamiltonian can be written in terms the occupation numbers of rotated lattice operators defined by where the function h(k) satisfies |h(k)| 2 = ε(k) 2 and reads It results in here we defined the energy bands by ε ± σ (k) = ±ε(k) with where a is the carbon-carbon distance (a ∼ 1.24Å). The energy resulting from a small variation on the occupation numbers on (28) at finite chemical potential µ reads where use have been made of the fact that the dispersion relation do not have spin indices to write the resulting expression in terms of variations of n ± k = n ± ↑k + n ± ↓k . To fix the ideas we consider only non-vanishing variations of the occupation numbers in the upper band, i.e. δn − k = 0, δn + k = 0. This can be done without loss of generality due to the symmetry of the system under the interchange of the positive and negative bands. We then have This has the form o the free term of (1), giving us one of the ingredients required by our step 1 defined above, namely the dispersion relation ε(k).
The band structure is shown in Fig 1a. The density of states is presented in Fig. 1b where a Van Hove singularity can be seen to be present for µ = ±t. The resulting FS's for different fillings are shown in Fig 1c.

Interactions in graphene
To complete the ingredients required by our step 1, we need the quasiparticle interaction function f (k, k ′ ). In what follows, for completeness and to set up our conventions, we briefly describe how to derive its expression to first order in a perturbative expansion [10,25].
We will consider density-density interactions, both onsite (with strength U ) and between nearest-neighbors (with strength V ), namely our interaction Hamiltonian reads where i, j stands for nearest neighbors, and the density operators n i = n i↑ + n i↓ refer to the original lattice operators a i , b i . One then computes the energy in the mean field approximation. The result can be written in terms of the mean field values of the occupation numbers of the rotated lattice operators that diagonalize H 0 , reading where N is the number of sites and (34) As advanced in the previous section, we will concentrate in variations of the occupation numbers that keep the total magnetization constant. In other words, we assume δn + ↑k = δn + ↓k . Similarly to the free part, the interactions between quasiparticles can then be written in terms of the variation of the total number n + k = n + ↑k +n + ↓k as The function f (k, k ′ ) is then obtained from the mean field value of the energy We have then completed step 1, getting the dispersion relation (29) and the interaction function (36).

Parametrization of the Fermi surface
Step 2' requires the parametrization of the FS, for which we need to study separately fillings lying above and below the Van Hove filling. In what follows we present the parameterized curves used throughout the calculations.

High energy sector: |µ| > t
We call high-energy sector the case in which t < |µ| < 3t, i.e. the region of fillings which lie above the Van Hove singularity. As can be seen in Fig.1c, for µ/t ∼ 3 the FS is approximately circular, while for values closer to 1 the FS takes a hexagonal form.
In this sector the FS can be parameterized as follows with ω H , φ H and the auxiliary function G(x) defined as Low energy sector: |µ| < t The low energy sector corresponds to fillings satisfying 0 < |µ| < t. In this case the FS consists of pockets centered at the vertices of the Brillouin zone, as can be seen in Fig.1c. By using the periodic identifications of the momentum plane, we see that only two of them are nonequivalent. In consequence one can describe the total FS as two pockets centered in the two Dirac points k ± = (0, ± 4π 3a ). For example, for the FS pocket centered in (0, 4π 3a ) we can choose a parametrization of the form and With the above parametrizations of the high and low energy sectors, we can compute the corresponding Jacobian evaluated on the FS according to (24), obtaining Following step 5 by means of an orthogonalization procedure, we obtain the new basis of mutually orthogonal functions {ξ n } n∈N .
In terms of this new basis and according to step 6, we compute the stability parameters χ i = ξ i , ξ i . These are functions on the space of parameters (µ, U, V ) that may become negative in some regions. If this is the case, in such regions the Fermi liquid is diagnosed to be unstable.
In our calculations, this last step was performed numerically, due to the complication of the integrals involved in the pseudo-norms {χ n } n∈N .

RESULTS
The phase diagram on the space spanned by the interaction strengths U , V and the chemical potential µ is shown in Fig 2. There we plot the instability region determined by the dominant unstable channels.
The method used in this work allows to explore all possible fillings and to draw a phase diagram for graphene systems valid both below and above half-filling. The advantage of our approach lies in the fact that it can be pursued systematically following the steps described in previous sections, studying separately each deformation mode of the FS.
For fillings around the Van Hove filling the Pomeranchuk instability is favored. For on-site and nearestneighbor Coulomb repulsion (U > 0 and V > 0) we find a region of parameter space where the system presents an instability. Near to the Van Hove filling our results are in agreement with those found using a mean field approach in [21].
On the other hand, we do not see any instability around half filling. This is in agreement with the results presented by Sarma et al in [22] for doped graphene, where the authors found that extrinsic graphene is a well defined Fermi liquid for low energies, within the Dirac fermion approximation. Moreover, this agrees with experimental results found using ARPES presented in [23], implying that graphene is a Fermi liquid for low dopings.
For attractive on-site interaction, the region where the instability is detected depends strongly of the nearest neighbors interaction strength.
We find that even the smoothest deformations of the FS, i.e. those described by lower modes in our orthogonal basis, present instabilities. Indeed, they cover most of the unstable region. In Fig.3 instability regions corresponding to the first modes are drawn. In Fig.3a the V = 0 plane is shown, the unstable mode corresponding to the colored region is the 0-th mode and the corresponding FS deformation is presented in Fig.3c. In Fig.3b and 3d the instability regions for the first modes are plotted in the U -V and µ-V plane respectively. In both Figures the χ 0 instability can be seen, and a new region appears where an instability of the χ 1 and χ 2 channels is present. This instability appears for positive values of the interaction strength and is closer to the Van Hove filling. This is in agreement with earlier results, where a Pomeranchuk instability at the Van Hove filling was found within a mean field approximation [21]. The interacting FS presented there has the same geometry than that of the corresponding deformation channel χ 2 shown in the left column of Fig.3c.
The results presented above confirm earlier predictions about the Fermi liquid behavior of graphene for the cases near Van Hove filling and near half filling. They also provide a more complete description of the phase diagram for the entire range of fillings in between these two limiting cases.

SUMMARY AND OUTLOOK
We have explored Pomeranchuk instabilities in graphene using a recently developed generalization of Pomeranchuk method. We have obtained the phase diagram of the dominant instability as a function of the on-site U and nearest neighbor V interaction strengths and the chemical potential µ (Fig.2). We analyzed several planes of the 3D phase diagram obtaining a good agreement with previous theoretical results and experimental findings (Fig.3).
The phase diagram makes apparent some interesting features of the system. For example no instability is detected at low energies. This behavior is noteworthy because in this sector a Dirac massless fermions approach can be used to describe the graphene layer in the absence of interactions. On the other hand, for energies close to the Van Hove energy, instabilities appear to cover a large region in the U -V plane.
The efficiency of the method shows up in the fact that the complete phase diagram is obtained by exploring a few number of modes. The introduction of higher modes does not enlarge substantially the region of instability but re-draw the details of its boundaries.
The method can be applied either analytically or numerically according to the complexity of the system under investigation. In the previously studied case [17], the analytical calculations were pursued up to the end, allowing us to draw the phase diagram. In the present case, the calculations were performed analytically up to the point of the evaluation of the instability parameters χ i , which involved complicated integrals that were then performed numerically. The method is suitable for direct application to numerical data encoding the form of the Fermi surfaces, like those obtained by application of ARPES.
The form of the method presented here is suitable for any two-dimensional lattice model at zero temperature. However, it does not consider instabilities arising from particles jumping between different disconnected pieces of the FS or from spin or color flipping. It can be easily extended to consider these effects, as well as to three dimensional systems, such as multilayer graphene, ruthenates, etc. The generalization to finite temperatures involves a different definition of the pseudo scalar product. All these extensions will be presented in a forthcoming work [27].

APPENDIX
Given the implicit curve defined by g(k) = 0, we can choose a parametrization in terms of a vector function k(t) = (k x (t), k y (t)) , −π < t < π , depending on an arbitrary parameter t, and defined so that ∀t : g(k(t)) = 0. In terms of this parametrization we want to prove that the Dirac δ function can be written as δ(g(k)) = dt |k(t)| |∇g(k(t))| δ (2) (k − k(t)) .
This is an identity in virtue of (48) if k y takes the place of x. Then we just proved formula (47) that we used during our calculation of the Jacobian of the change of variables evaluated in the FS.