Symmetry breaking and lattice kirigami: finite temperature effects

Recent work has analysed how deformations due to the insertion of a defect in a flat hexagonal lattice affect the ground state structure of an interacting fermion field theory. Such modifications result in an increase of the order parameter in the vicinity of the defect and can be explained by a kirigami effect, that is the combined effect of the curvature, locally introduced by the deformation in the lattice tessellation, and of a synthetic gauge field induced by the boundary conditions along the cut, performed to introduce the defect. In this work, we extend the formalism and previous results to include finite temperature effects.


I. INTRODUCTION
Quasiparticles in quantum materials are influenced by the configuration of the crystals in which they move. In realistic crystals, in fact, the unavoidable presence of defects induces an effective change in the topology and the geometry of the lattice, with drastic consequences on quasiparticles propagation. The rapid development of this new subject in the context of quantum fields in condensed matter systems has sparked the interest in a rejuvenation of methods and results of semiclassical gravity. In the continuum limit, one can think of the lattice in which these particles move as a curved background, and apply some of the well-established techniques of quantum field theory in curved space [1,2].
Most of all, graphene [3] is the example that recently has attracted attention on both the theoretical and applicative grounds. Graphene is the archetype of the quantum material realising an emergent system of relativistic fermions moving in a bidimensional space. Since its first synthesis in 2004 [4,5] (but even before, following some pioneering experiment in the 60s [6][7][8]), much work has been devoted to the methodic engineering of those conditions that alter the conductivity properties of graphene. Fermions self-interactions [9,10], impurities arising in the fabrication process [11,12], the substrate on which the monolayer graphene sample lies [13], all of them could in principle be responsible for the appearance of a mass gap in the energy spectrum of the quasiparticles propagating on the carbon sheet. The case of crystal defects is particularly interesting: several kinds of structural defects have been isolated and investigated since early experiments, among the others Stone-Wales defects (a pair of pentagons, separated by a pair of heptagons) [14], mitosis defects (a pair of heptagons, separated by a pair of pentagons) [15], single and multiple vacancies, line defects (chains of pentagons, heptagons and hexagons gluing two patches of the crystal) [16], or out of plane carbon adatoms. Quite surprisingly, the competition between the topological features and curvature arising after * flachi@keio.jp † vincenzo.vitagliano@keio.jp the introduction of the defects has shown definitely nontrivial aspects [17]. Similar considerations have been raised for more complex crystalline structures: for example Weyl semimetals [18][19][20][21][22][23][24][25][26][27], three-dimensional topological materials whose low energy excitations are Weyl particles, that is massless Dirac particles with fixed chirality. Novel developments have enlightened the response of these materials to geometrical deformation and/or electromagnetic stimulation. In such a context, classical Einstein-Cartan geometry (see [28] and reference therein) provides the best framework within which modelling the effective theories that describe exotic quantum phases, with torsioninduced defects contributing in a critical way to current anomalies [29][30][31][32][33][34][35].
Fermion conductivity of quantum materials is in general sensible not only to the intrinsic conformational properties of the hosting lattices but also to the temperature of the system and eventually to the occurrence of inhomogeneous and anisotropic phases. Electrical conductivity in graphene has been the object of investigation in the framework of the Dirac model, at arbitrary values of temperature and chemical potential, in both the gapless and gapped configurations (see, e.g., [36]). Here we are interested in studying the interplay between the effects due to finite temperature and finite density with those related to geometrical issues on a system of interacting particles lying on a hexagonal (graphene-like) lattice.
As previously mentioned, four-fermion interactions are amongst the possible concause for the emergence of a non-vanishing mass gap in graphene: collective excitations arise as a consequence of some symmetry breaking mechanism (for graphene this could happen for example for breaking both Z 2 and the discrete sublattice symmetry, which leads to staggered magnetisation). However, the occurrence of symmetry breaking in a nonflat manifold is modified by the presence of curvature [37], with the Ricci scalar entering the gap expression in terms of an effective mass, thus contributing to the total mass of the condensate. The problem that this paper will address can be formulated as how the gap shift determined by geometry affects the phase diagram of the interacting field theory on the honeycomb lattice. Previous results on QCD phase transitions in a strong gravity environment [38,39] have tracked the route for the development of a formalism based on the use of the effective action and zeta function regularisation to study particle condensation under the influence of external fields in a nonperturbative way. In the next sections, we will see how this proposal can be extended to the case of symmetry breaking on engineered curved lattices.

II. THE MODEL
A simple, while still effective, description of strongly correlated particles on a lattice is given by the singleband Hubbard model [40]: being u † , u and n respectively creation, annihilation and number operators, the Hubbard Hamiltonian reads the first term expresses the kinetic energy of the system: a particle with spin projection σ (= ↑, ↓) is destroyed in the l-site and created in the j-site; in the Hubbard approximation, this particles' hopping is allowed only between two adjacent sites (this is the meaning of the symbol < j, l >), corresponding to the positions of two first neighbour atoms on the hexagonal lattice; this assumption is justified by the exponential drop-off of the particles' wavefunctions profiles, with the hopping energy scale t set by the overlap of two contiguous wavefunctions. The second term in (1) controls the Coulombian repulsive interaction energy between two particles with opposite spin occupying the same site j, with a coupling strength set by the parameter U > 0. Finally, the third term, introduced by a non-vanishing chemical potential µ, describes the particles filling of the lattice sites. The Hubbard Hamiltonian is invariant under a global SU(2) spin transformation and a U(1) (charge) redefinition of the one-particle wavefunction. On top of these, further symmetries arise if the space where the particles live is a lattice with a bipartite structure. A bipartite lattice is obtained from the union of two interpenetrating sublattices A and B (e.g. the red and blue dots triangular sublattices of Fig. 1) such that the first neighbours of an A site are all B atoms. Hexagonal and square lattices are bipartite lattices. Triangular lattices are not. Bipartite lattices' energy is trivially minimised by Néel states, where the spins of one sublattice are all parallel among them and antiparallel to the spins of the second sublattice. In lattices which are not bipartite, instead, the Néel state is generally degenerate with other classical configuration obtained from the former by a local spin flip; these systems are said to be geometrically frustrated: half-filled lattice configurations, i.e. systems for which each site in the lattice has one and only one fermion, always have at least a pair of contiguous sites occupied by particles spinning in the same direction.
When the background on which the particles move is a bipartite lattice, the Hubbard Hamiltonian enjoys additional symmetries: it is invariant under a particle-hole transformation, that is under the change of creation into annihilation operators (and vice versa) Moreover, its spectrum remains unchanged under a sign flip of the hopping parameter t.
The influence of temperature and geometry on the occurrence of symmetry breaking is more conveniently studied in a covariant functional-integral formulation of the Hubbard model associated partition function, which can be straightforwardly achieved [41][42][43] bosonising the density-density product in the interaction term of (1). An explicitly SU(2)-invariant Hubbard-Stratonovich transformation of the interaction term may be obtained, for spin-1/2 fermions, by introducing an arbitrary unit vector n along the spin-quantization axis of the particles, such that where Σ = (σ x , σ y , σ z ) is the usual Pauli vector. Introducing two auxiliary fields φ and ∆, the density-density operator can be rewritten as The effective Lagrangian for the fluctuations is then derived integrating out the fermionic degrees of freedom.
FIG. 1. The hexagonal honeycomb lattice with the two triangular sublattices (blue and red spheres).
The procedure here described is entirely general; however, to ease the calculation, we will consider only one scalar order parameter, φ j , breaking at the same time the Z 2 and the discrete sub-lattice symmetry of the halffilled Hubbard model on a hexagonal lattice. This order parameter can be identified with the staggered magnetisation -the conjugated of the density operator -defined as the net difference between the number of negatively and positively spinning fermions, φ j ≡ n j↑ − n j↓ .
For sake of compactness, it is possible to introduce a new wave function ψ σ such that ψ T , where X = A and B refer to the two sublattices, # = 1 and 2 are the two inequivalent Fermi points K # of the half-filled honeycomb lattice, and ψ X,# σ are the inverse Fourier transforms of the sublattice annihilation operators u X,σ (p + K # ). With these assumptions the effective Lagrangian finally reads expressed in terms of the usual flat space γ-matrices, γ a ; here λ is a constant, proportional to the coupling strength parameter U up to an irrelevant factor, and α ↑,↓ = ±1.

III. KIRIGAMI AND LATTICE GEOMETRY
The Japanese term kirigami describes a particular origami technique which, apart for folding the paper, allows for cutting parts of it. Two-dimensional carbon allotropes have recently been shown to be perfect candidates for engineering robust three-dimensional microscale structures with tunable mechanical properties (see for example [44]), making of them a sort of "lattice kirigami". The most straightforward procedure to obtain such kirigami is typically via defect insertion in the (honeycomb, in the present case) lattice, as shown in Fig.  2. Starting from a locally flat lattice where a hexagonal cell is marked as the centre, one can design two kinds of defects, bringing to two different geometries: reducing the number of the central cell sides of N > 0 sides corresponds to cut and pulling out a N π/3 section of the lattice; gluing together the two sides results then in a conical configuration with positive curvature. On the other side, augmenting the number of sides (which in our notation corresponds to consider a negative N ) of the central cell is equivalent to pull into the lattice a |N |π/3 sector, resulting in a saddle geometry with negative curvature. We will only consider lattices with an even number of defects, preventing possible frustration effects from veiling or even obscuring the main quantum mechanisms here studied.
The kirigami procedure, as stated above, alters the topological and geometrical features of the material sheet: the very same fermions propagation is sensitive to these conformational properties of the underlying lattice. As a result, the study of the vacuum structure and the symmetry breaking behaviour exhibited by a continuous field theory like (4) defined on the curved lattice, relies on the description of an effective continuous geometry of the lattice itself.
In the continuum limit, this configuration of the lattice can be described by a conical metric [45][46][47]. Away from the central defect (corresponding to a conical singularity), the associated spacetime is accurately described by a locally flat geometry, ds 2 con = dt 2 − dr 2 − r 2 d θ 2 , but with an important caveat: it is not globally Euclidean, since the angular coordinate does not run on a 2π circle; instead, 0 ≤ θ < 2π − ∆, with ∆ > 0 (∆ < 0) being the deficit (excess) angle: surfaces at constant t are cones (saddles), not planes. It is worth to stress that in the case under examination (lattice with a hexagonal base) ∆ does not take continuous values: it is clear that it should be ∆ = N π/3, with |N | even to avoid the frustration of the system mentioned above. Note also that there is a constitutive bound on N : while the only positive value allowed is N = 2 (corresponding to the square defect), negative values are allowed up to |N | ∼ 16, on top of which the geometrical structure is supposed to collapse onto a helical conformation [48].
The conical metric, ds 2 con , is singular on the tip, where it is not possible to introduce a tangent space and consequently, is not possible to calculate the curvature. This singularity is only an artefact of our mathematical description of the cone: the deformed lattice does not present any singularity in correspondence of the defect. A better approximation of the physical lattice is achieved smoothing the geometry into that of a snub-nosed cone (resp. saddle), with metric ds 2 regulated by a parameter . The biggest is the regulator , the smoothest the tip of the cone. The proper conical metric ds 2 con is recovered as the limiting case of a sequence of such smoothed cones, when ds 2 con = lim →0 ds 2 (see Fig. 3). In the lattice language, the regulator is never vanishing, rather it parametrises the minimum physical distance between adjacent atoms.
Let us introduce a new angular coordinate θ = θ/α, with α = 1 − ∆/2π, such that 0 ≤ θ < 2π. The line element of the regularised tip can then be written as where the smooth function f (r, ) satisfies the following asymptotics: i) lim r→0 f (r, ) = α 2 and ii) lim r f (r, ) = 1. The simplest choice implementing these requirements is the sequence of spaces with metric The covariant extension of the model (4) to the curved spaces (6) is obtained under the minimal coupling ansatz : the flat metric is substituted by the metric tensor g µν defined in (6) and such that g µν = e µ a e ν b η ab , being e µ a the vielbein and η ab the Minkowski metric; the ordinary partial derivative is replaced by the covariant derivative and the flat γ-matrices by γ µ = γ a e µ a , which are the γ-matrices in curved space.
There is still one geometrical issue arising in the kirigami procedure. The cut-and-paste technique to induce lattice defects changes the topology of the system. One can show [49] that the boundary conditions satisfied by the spinors along the cut read as what is evocatively called the Moebius stripe conditions: rotation around the defect brings to a configuration where the two triangular sublattices are exchanged, while one more rotation brings the system back to the initial state, where R = i 0 σ 2 −σ 2 0 = −γ 5 (last equality holds in the standard planar representation, where γ 0 is diagonal, but changes in other representations of the Clifford algebra). Note that R anticommutes with all the γ-matrices. Using a gauge-like transformation mapping ψ into a new field ψ such that it is possible to re-arrange the Moebius boundary condition (7) to resemble the flat space (antiperiodic) condition ψ(r, θ + 2π) = − ψ(r, θ) .
The effective Lagrangian (4), under the fields redefinition (8), acquires an extra contribution in the form of a (constant) effective gauge connection term 1 , B µ ≡ −δ θ µ N 4 R, so that now the bosonised Hubbard model reads where D µ = ∇ µ + ıB µ . In what follows, we will drop all the tilde's from the previous equation to shorten the notation, having in mind the caveat that spinor fields have already undergone a Moebius stripe rotation.
In order to understand the phase diagram of the lattice kirigami we need to study the behaviour at different temperature of the staggered magnetisation parameter φ, which can be calculated solving numerically the effective equations of motion for the system. Keeping this strategy in mind, we start by evaluating the Euclidean effective action for the model (10) on the Wick-rotated curved background (6), As it is, this expression of the effective action does not allow to obtain directly the curvature counterterms to the one loop divergences of the Dirac operator (the second term in (12)); this is due to the circumstance of O ≡ (ıγ µ D µ + α σ φ + µγ 0 ) being a linear differential operator of the first order, rather than being of second order. A possible way out is to apply standard techniques and square the Dirac operator: the eigenvalues of the squared (Dirac) functional determinant on the Riemannian spin manifold are then related to those of the sum of the spinor Laplacian with the space curvature, through the Schrödinger-Lichnerowicz-Weitzenböck formula, / D 2 = ∇ * ∇ + I · R/4 (see for example [57]). Following this prescription, the one loop effective action (12) is rewritten as where it should be taken into account that the d'Alembertian operator,ˆ , is built out of the total covariant derivative including the contribution from the effective gauge connection and calculated on the Euclidean version of the regularised metric, ds 2 (see also [38,58] for further details); in the last line we have taken advantage of the mode expansion (11) in the Matsubara frequencies. The expression of the effective action (13) can be obtained through different techniques, usually always quite cumbersome. In the present case, we will use zeta-function regularization (see for example [1,[59][60][61][62][63]). Given an operator O, its functional determinant is the product of the N eigenvalues κ N , such that where l is some renormalization length. Defining the function it is possible to rewrite the expression (14) as provided that ζ and its first derivative are regular in zero. In three dimensions, ζ(s) converges for [s] > 3/2, and can be analytically continued to a meromorphic function of s for [s] < 3/2; in particular it will be regular in zero, and will encounter simple poles for s = (3/2 − p), p ∈ N. Under these assumptions, and using the integral representation of the function Γ(s), it is possible to relate the zeta-function to the Mellin transform of the heat trace, which in the case of the differential operator in (13) reads In curved space, the heat kernel of a second order differential operator like the one appearing in (17) has a particularly enjoyable proper time expansion, exhibiting a re-summation, among the curvature invariants contributions, of all the power terms in the scalar curvature R [64][65][66][67], where the traces are taken on the Dirac indices and the first three D (k) σ coefficients read with In order to calculate the zeta-function (17) we will follow the strategy delineated in [38]: the sum on the quantum number n can be expressed in a more convenient way observing that (see also [68] (−1) n e − β 2 n 2 4τ cosh(βµn) ; (20) using (20) and (18) in (17) one finds, after some lengthy calculation, an expression for the zeta-function. In odd dimensions, the analytic continuation of ζ(s) to s = 0 vanishes, while the first derivative turns out to be where K's are the modified Bessel functions of the second kind (which have the nice property to be exponentially suppressed) and being the digamma function (the logarithmic derivative of the gamma function Γ). Considering only terms up to the second order in the heat kernel expansion, (21) reads which, given that ζ(0) vanishes, straightforwardly gives the logarithm of the functional determinant, and hence equal to the effective action. Substituting into (14), one can finally use the obtained expression to calculate the effective action corresponding to the theory with Lagrangian (4) and metric (6). The effective action can now be numerically evaluated. As an example, Figs. 4 and 5 show the results for zero chemical potential, positive curvature (cone) at different values of the temperature. The two figures in Fig. 4 show, respectively, the trends of the effective potential close to the apex of the regularised cone (left panel) and far away from the defect (right panel), as a function of the temperature. Note that, as expected, increasing the temperature corresponds to pushing the system toward a disordered phase. The numerical integration of the effective equations of motion finally gives an explicit profile for the order parameter, Fig. 5: the approach to the central defect enhances the symmetry breaking and the formation of the condensate, further confirmation of the effect discovered in [17], with temperature favouring the phase transition.

V. CONCLUSION
In quantum field theory, the Coleman-Weinberg mechanism [69] explains how symmetries may spontaneously break as a result of quantum effects, predicting a phase transition (of the first order in scalar electrodynamics) from a broken to a restored symmetry phase as the mass is increased. Once the background geometry underlying the field theory is curved, the same mechanism leads to expecting a similar transition when the scalar curvature is increased. However, this is not the whole story, and when the topology of the background becomes nontrivial, new effects may take place. An example of the sort has been studied in Ref. [17], with the background geometry emerging from the continuum limit of a deformed hexagonal lattice, with the deformation being due to the presence of a defect in the lattice tessellation. This "kirigami-engineered" background is quite special, since aside for altering the curvature it also modifies the boundary conditions that fields should obey when circulating the defect. These non-trivial boundary conditions can be mimicked by a synthetic gauge field localised near the apex and compete with the curvature in altering the ground state, inducing condensation close to the defect. A relevant question is how additional external conditions change the picture. To address this problem, we have adapted to the present case the imaginary time formalism and zeta function regularisation techniques to arrive at a partially re-summed form of the effective action, whose minimisation has been carried out by numerical approximation. We have studied what happened when effects of finite temperature are included, resulting in a technically non-trivial modulation of the order parameter, as illustrated in Fig. 5. Given a certain temperature, the condensate is always enhanced when approaching the defect, as a consequence of the kirigami effect.
An intriguing possibility is to tune the temperature to obtain a more complicated pattern of the condensate as a function of the distance from the defect. The black curve of Fig. 5, for example, shows the possible profile of a condensate obtained decreasing the temperature as stepping away from the defect. Here, the interesting aspect is in the value of the condensate, that increase with dis-tance, rather than decreasing. The situation pictured in this way, however, is valid only as a first approximation. Temperature anisotropies, in fact, induce extra variation in the local density of the condensate and could make the effect of inhomogeneities (and hence of the contribution of the chemical potential part to the effective action) non-negligible.