Tunning the tilt of a Dirac cone by atomic manipulations: application to 8Pmmn borophene

We decipher the microscopic mechanism of the formation of tilt in the two-dimensional Dirac cone of $8Pmmn$ borphene. In our ab-initio calculations, we identify relevant low-energy degrees of freedom on the $8Pmmn$ lattice and find that these atomic orbitals reside on an effective honeycomb lattice (inner sites), while the high-energy degrees of freedom reside on the rest of the $8Pmmn$ lattice (ridge sites). Integrating out the high-energy atomic orbitals, gives rise to remarkably large effective further neighbor hoppings on the coarse grained"honeycomb graph"of inner sites that determine the location and tilt of the Dirac cone. Molecular orbitals -- that can be modified by atomic manipulations -- are responsible for the creation of further neighbor edges on the honeycomb graph that controls the tilt of the resulting Dirac cone. This leads to an effective tight-binding model on a parent honeycomb graph that facilitates numerical modeling of various effects such as disorder/interactions/symmetry-breaking for tilted Dirac cone fermions. Since the tilt is a proxy to spacetime metric, our result offers a robust perspective on the fabrication of desired emergent spacetime structure and synthesis of geometric forces at ambient conditions that are likely to enhance our control on the movement of electrons in electric/optical devices.


I. INTRODUCTION
The marriage between the mathematical concept of "topology" and quantum materials lead to the birth of "topological materials" and developments that followed afterwards. Can the (local) "geometry" play a similar role in combination with quantum materials? In the same way that the appearance of Dirac fermions in graphene and other Dirac materials can be attributed to an emergent Minkowski spacetime structure at distances much larger than the atomic distances, quantum materials with "tilted Dirac cone" can be naturally associated with an emergent spacetime (not merely the space) metric. In these class of materials the tilt is a proxy to spacetime geometry. We use Geometric Quantum Material (GQM) to refer to them.
Every periodic structure in quantum condensed matter is mounted on a mathematical object called lattice. Irrespective of which atoms one wishes to place on the sites of a given lattice, they come in 230 possible structures [1]. The presence of lattice breaks translation and rotation, and some times parity and/or time reversal invariance of the vacuum [2]. The lattice therefore breaks the Poincaré group [3] thereby the connection between spin and statistics of the particles is lost and therefore one may have fermions with integer spin that is enforced by the irreducible representations of its space group (SG) [4].
Is there any other interesting consequence that can be associated with the underlying SG? To set the stage for answering this question, let us start by a simple, but profound observation [2]: Consider a simple quantum mechanical hopping process on a lattice and think about the wave equation in the continuum limit. On the square lattice the Hamiltonian becomes −h 2 ∇ 2 /(2m * ), while on the honeycomb lattice it becomes ihv F σ.∇, namely the massless Dirac Hamiltonian [5]. Therefore despite that in the continuum limit, the lattice spacing is immaterial, but still remnants of the microscopic symmetries of the pertinent SG are manifested in the long-distance behavior and decide whether the structure of the ensuing spacetime is Galilean or Minkowski. This can be regarded as an example of metamorphosis of atomic scale symmetry (SG) to a long-distance geometrical structure. The duality between atomic scale SG symmetry and long distance geometry is a central to understand the geometric structure in GQMs.
What is the microscopic mechanism of such a duality between the microscopic structure and the long-distance (low-energy) characteristics in more complicated lattice structures? The long-distance behavior of interest for us is the tilted Dirac cone band whose Hamiltonian is ∝ (σ.∇) + σ 0 ζ.∇ [6-8, 10, 11, 19] where σ 0 is the 2 × 2 unit matrix. A vector-looking quantity ζ here determines the tilt of the Dirac cone. But at a much deeper level, it can be encoded into an emergent spacetime metric ds 2 = −(v F dt) 2 + (dx − ζv F dt) 2 [10,[12][13][14][15][16][17]19], where v F is the Fermi velocity. Therefore the tilt of a Dirac cone is actually a proxy for an emergent solid-state spacetime structure. Since the density of states is enhanced by a 1/ 1 − ζ 2 , the tilted Dirac cone materials can generically give stronger responses to external stimuli 1 . Therefore it is crucial to understand the atomistic mechanism of the formation of the tilt in order to utilize it. In this work we decipher the microscopic mechanism that determines the tilt parameter ζ. For this purpose we focus on the SG number 59 that describes the 8P mmn structure 2 of borophene [19][20][21]. But the logic employed here is generic and is applicable to other SGs, and can even be used to discover new examples of GQM. We show that in this structure, the lowenergy and high-energy degrees of freedom are nicely separated into two sublattices denoted by gray and teal circles in Fig. 1(a), respectively. When the high energy sites available in this particular SG are integrated out, the underlying honeycomb lattice will be promoted to a "honeycomb graph". Molecular orbitals play significant role in attaching new "graph edges" to the simple underlying honeycomb lattice. This graph has the ability to lead to a controlable tilt and hence spacetime structure in the long distance. On such a graph, even if instead of fermions one places a circuit, again a tilted Dirac cone can be obtained [22]. Therefore GQMs themselves can be emulated by circuits with a larger degree of tunability.

A. Protection of the Dirac node
The relevant orbitals in the boron (as well as C) atom are 2p orbitals. The possibility of formation of sp 2 and sp 3 hybridization establishes the honeycomb lattice [23], and structures such as 8P mmn that involve buckled honeycomb networks as natural lattices for these atoms. As can be seen in Fig. 1(a), there is a backbone (buckled) honeycomb sub-lattice denoted with gray circles that are called inner (I) sites. The rest of the lattice sites are called ridge (R) sites and are denoted by teal circles. First principle calculations indicate that the resulting Dirac cone is tilted [24,25] which is shown as red (low-energy) band in Fig. 1(b) and the three dimensional reconstruction of the band structure is shown in panel (c). The first thing that the 8P mmn SG implies   [26]) where the two-dimensional irreducible representations are protected by non-symmorphic elements [26]. Then the compatibility relations gives the qualitative band picture of cat's cradle shape [27] shown in panel (d).
For pedagogical details of the derivation of this figure with group theory methods see the supplementary information (SI). Of course the protection of the tilted Dirac cone by the underlying SG is interesting, but it is not essential for the main purpose of this paper, namely the "tilt" of the Dirac cone. The symmetry considerations do not explain why and how the tilt parameter ζ is formed. In order to "manipulate" the tilt of the Dirac cone at will, one needs a microscopically detailed understanding of the root cause of formation of tilt in the Dirac cone. To achieve this, we need to identify the low-energy degrees of freedom that give rise to the tilted Dirac cone dispersion. For this purpose in Fig. 1(e) we have plotted an orbital projected representation of the band structure that resolve the contribution of atoms 2 (I) and 5 (R). As can be seen, the dominant contributions to the tilted Dirac dispersion comes from the p z orbitals of the I atoms (top left plot in panel (e)). The R atoms also contribute via their p z and p x orbitals (bottom left and center). The remaining two boxes on the top row indicate little contributions from the p x and p y orbitals of the I atoms that arises from the buckling of the honeycomb sublattice.
In panel (a) of Fig. 2 3: (a) Virtual hopping via the ridge site 7 generates an effective hopping t p between 2 and its third neighbor site 3 -i.e. site 3 in the adjacent unit cell in Fig. 2 -(b) The odd parity combination of 2 and 3 remains decoupled, but the even parity combination hybridizes with ridge site to gain energy. hoppings through R sites. There is only first neighbor direct hopping between the I sites denoted by black arrows. The appropriate atomic orbitals involved in forming the above microscopic t ij hoppings can be extracted from the DFT calculation. Working in a gauge that all t ij 's are real, the hermiticity implies t ij = t ji . Panel (b) depicts the BZ of original 8P mmn lattice. The insight from the projected bands in the left and center columns in the second row of Fig. 1(e) is that both p x and p z orbitals of the R atoms are of comparable importance and must be incorporated into the atomic scale computation of the t 65 (dark green) and t 87 (light green). Similarly the first column of Fig. 1(e) suggests that the p z orbitals of R and I atoms dominate in t 36 = t 35 hopping process. Tab. I shows the calculated values of these hopping parameters for pristine borophene (B 8 ) and C-doped borophene B 6 C 2 using Wannier function [16,17]. The purpose of substituting C for B is to study its effect in the tilt of the Dirac cone. The substituted carbon dimers are placed in R and I positions, respectively. Using these parameters one can re-construct the ab initio bands with four p z orbitals of four inner sites, and eight p x and p z orbitals of ridge sites. The nice coincidence of the original DFT bands with Wannier-interpolated bands in Fig. S9 (see SI) shows that the obtained t ij values are reliable and the used atomic orbitals are adequate.
Although such a atomic picture might be satisfactory if one wishes to focus on the low-energy features of the tilted Dirac cone around the Fermi surface, but still working with a 12-band Hamiltonian is neither convenient, nor the essential long-range physics depends on so many short-distance details. To achieve an effective two-band model, we need to decimate the 8P mmn lattice into an effective honeycomb graph shown in Fig. 2(c) where two possible ways of representing its BZ are depicted in panel (d). On such a coarse-grained lattice, the virtual atomic hopping paths will be replaced by effective hoppings of the same color in panel (c). The connection between these effective hoppings and atomic hoppings t ij is a nice example of renormalization that will be discussed now.

B. Renormalization via molecular orbitals:
Anderson in his book maintains that renormalization is one of the pillars of condensed matter physics [31]. In this section we will show that the same concept is encoded into a local quantum chemistry of the 8P mmn borophene in a remarkable way. Since we are not interested in higher energy features that take place away from the tilted Dirac node, we build an effective picture based on Fig. 2(c) where the "effective" hoppings must be evaluated from the atomic scale data in the upper part of Tab. I. The low-energy degrees of freedom are dominated by the p z orbitals of the I sites. So we decimate the 8P mmn lattice by elimination of the R sites. In doing so, the effective coarse-grained system becomes the honeycomb graph in Fig. 2(c). The nearest neighbor hoppings denoted by black arrows take place within the low-energy subspace and are responsible for the formation of a parent Dirac cone from the p z orbitals of the I sites. Slight anisotropy in the black arrows is known to shift the location of the Dirac cone [32] within the rhombus BZ of panel (d) that we ignore here. Green/red hopping processes are associated with the second/third neighbors of the honeycomb backbone that originate from the corresponding atomic process t ij of panel (a) via the renormalization.
Let us see how the atomic processes t ij in Fig. 2(a) are related to the effective hopping processes in Fig. 2(c). For example consider the simplest third neighbor R sites 2 and 3 of adjacent unit cells in Fig. 2(a) that becomes possible via atomic hoppings t 27 and t 37 of Fig. 2(a). The hopping process via the R site 7 is depicted in Fig. 3(a). Assuming an on-site energy offset 2∆ for the R sites with respect to I site, an electron starting at the site 2 virtually hops to R site 7 and then returns to low-energy sector at site 3. Through this process it gains the following energy  Even in the limit of ∆ → 0 this gives an energy lowering − √ 2|t 27 | that can be regarded as effective hopping between the third neighbor sites 2 and 3 of the I sublattice. To intuitively understand this formula, note that in the absence of the I site, t 27 = t 37 = 0, and hence both even and odd combinations of third neighbor atomic orbitals |2 ± |3 remain inert. But once the inner site is present, the atomic hopping t 27 causes a coupling of |7 with the even-parity state |2 + |3 , leaving the odd parity state |2 − |3 decoupled at zero energy. The coupling with the even parity combination provides a channel to lower the energy given by the effective hopping in Eq. (1). Furthermore, in the limit of ∆ |t 27 | the above formula reduces to the perturbatively appealing form −2t 2 27 /∆. For pedagogical details of the above computation, please refer to section II of SI. As a second example, consider the I site 3 in Fig. 2(a) and the same site in the unit cell above it (that we denote by 3 ). These sites will be second neighbor on the coarse grained lattice and there are two hopping pathways 3 → 6 → 5 → 3 and 3 → 8 → 7 → 3 connecting them that are denoted by light and dark green dashed lines in Fig. 2(a). The process of the calculation of the effective hopping amplitude 3 → 3 is similar for the above two paths. So we focus on the first one. In principle one must consider the energy gained by the lowest molecular orbitals formed by the above chain of sites. This has been done in section III of SI. But the end result allows for a nice interpretation depicted in Fig. 4: First, due to hopping t 65 , the R sites 5 and 6 form bonding and anti-bonding molecular orbitals denoted by |φ 1 and |φ 3 in Fig. 4(a). Then as depicted in panel (b), an electron gains energy by virtually hopping via the anti-bonding orbital |φ 3 whose energy is now offset by t 65 + 2∆. The hoppings connecting site 3 and 3 to |φ 3 are t 65 / √ 2 where the √ 2 factor comes from the normalization |φ 3 = (|φ 5 −φ 6 )/ √ 2. In the final step as shown in Fig. 4(c), the even-parity combination of |3 and |3 is mixed with |φ 3 to give an energy gain E even − = (t 56 /2 + ∆) − (t 56 /2 + ∆) 2 + t 2 63 . A similar contribution arises from the second path. Adding the two contributions we obtain the total renormalized hoppingt of Fig. 2(c) as In the above formula, the hopping parameters t 65 and t 87 are dominantly contributed by the p x and p z orbitals as the intensity of the p x , p z orbitals in the second row of Fig. 1(e) is dominant, while p y is faint. Similarly the two other renormalized hopping parameters t x andt can be computed. For details please refer to SI where we have shown that the elimination of higher energy states actually corresponds to the above simple molecular orbital analysis. This is how the renormalized parameters of the coarse grained honeycomb graph model in the bottom part of Tab. I are computed from the ab initio data of the top part of the table. Note that if instead of 8P mmn structure, we had a simple honeycomb lattice (such as in graphene), such a large values of second or third neighbor hopping given in Tab I would be unthinkable as hopping between the atomic orbitals exponentially decays with distance. Therefore the virtual hopping via the ridge sites attaches a great importance to them as providers of channels for energy gain and ultimate formation of renormalized hoppings on a coarse grained lattice of inner sites. Furthermore, this is an example of how the molecular orbital play a significant role in the formation of longer range hoppings on the decimated 8P mmn lattice.

C. Effective coarse-grained model
Now that we have identified the p z orbitals residing at the inner sites as the low-energy degrees of freedom, and have computed various renormalized hoppings between second and third neighbors, we are ready to write down a physically clear low-energy effective model that will straightforwardly demonstrate how the renormalized parameters t p , t x ,t andt can provide information about the position and the amount of the tilt. The first neighbor hoppings denoted by black arrows in Fig. 2(c) and are present even when there are no R sites. In the twodimensional Hilbert space of A and B sublattice of such a honeycomb lattice, these hoppings contribute the usual off-diagonal term F 0 (k) = − α=1,2,3 e ik.δα to the Bloch Hamiltonian, where the sum over α runs over three first neighbors in Fig. 2(c). This is responsible for the formation of a pair of upright Dirac cones similar to graphene. Slight anisotropy of the honeycomb lattice of R sites will amount to a shift in the location of the Dirac cone [32] which is irrelevant for our purposes. For concrete calculations, we assume a regular effective honeycomb lattice of bond length a for which the Dirac nodes in the rhombus BZ are at K ± b = 2π 3a (î ± 1 √ 3ĵ ) [23]. Next let us consider the third neighbor hoppings denoted by red arrows in Fig. 2(c) labeled by t p (p for pseudo, as they give rise to pseudogauge fields that shift the location of the Dirac node) and t x (x to emphasize the role of p x orbitals of R sites). They contribute another off-diagonal term denoted by = F xp (k) = t p e 2iakx + 2t x e −iakx cos √ 3ak y to the Bloch Hamiltonian. Expanding this form factor around the Dirac nodes K ± b above, gives rise to (i) a shift ∆k D = 2(t p − t x )/(∓3at ± at x ) of the Dirac node, where t is the nearest neighbor hopping (assumed to be 1) and (ii) anisotropy in the Fermi velocity given by 3t . Therefore the first and third neighbor hoppings establish the location of the (still upright) Dirac cone and determine its Fermi velocities. Now let us focus on the second neighbors (green arrows) in Fig. 2(c) that are the root cause of the tilt formation [7]. Since these hoppings are driven via virtual hopping through different arrangements of molecular orbitals, there are two types of them denoted by solid (t) and dashed (t) green lines. These hoppings being second neighbor, connect two atoms on the same sublattice, and therefore contribute to the diagonal terms, namely AA and BB components of the effective 2 × 2 Hamiltonian matrix. Among the matrices σ µ with µ = 0 . . . 3, only σ 0 and σ z can contribute diagonal terms. So now one has to decide whether these diagonal terms come with the same sign (σ 0 term → tilt) or opposite signs (σ z term → gap). There are two ways to see that this term must be proportional to the σ 0 : (i) Analysis of the irreducible representations and compatibility relations of the original 8P mmn structure in Fig. 1(d) shows that the crossing of the red bands is protected by the glide elements of the 8P mmn lattice (see SI for details). Since the effective theory has to obey this protection against gap opening, the σ z term is ruled out. (ii) Consider the renormalized lattice itself and focus on the solid green line in Fig. 2(c). If the hopping between 1 and 3 in Fig. 1(a) contributes to AA term, the hopping between 2 and 4 contribute to BB term. Both these contributions arise from the p z orbitals of these atoms via intermediate hopping through p z orbital of atom 5. Apparently for p z orbitals the "northwest" (1 → 5 → 3) and "northeast" (4 → 5 → 2) hoppings are identical. Similar arguments holds for the dashed green line in Fig. 2(c). In this case AA (BB) term is generated by 3 → 6 → 5 → 3 (2 → 6 → 5 → 2) path. Again the p x orbitals of the 5, 6 atoms symmetrically connect 3 → 3 and 2 → 2 (remember 3 is the same as 3 but in adjacent unit cell. Similarly for 2 .), thereby giving identical AA and BB terms in the effective Hamiltonian.
Therefore the effective Hamiltonian becomes Taylor expanding the diagonal f tilt term around the Dirac node formed by off-diagonal terms gives, where ± corresponds to the valley around which the expansion is performed.
The above equation indicates that the tilt arises from the difference of the second neighbor hoppingst andt that in turn are generated via the p x and p z orbitals of the R atoms. An immediate suggestion of the above model is to (partially) replace the R-site boron atoms by carbon atoms to see whether the tilt is changed or not. In Fig. 5 we have replaced two of the boron atoms with C atoms. For this purpose there are two choices: (i) To place a carbon dimer on the R sites as in the panel (a) or (ii) to place the carbon atoms in the I sites as in panel (f). The results are summarized in Tab. I. For case (i), the location of the Dirac node is shifted towards the Γ pint and its tilt increases. Placing carbon atoms in the R sites shifts the 5 and 6 sites to higher energies, thereby generating largert andt (see Tab. I) that ultimately increases the tilt from the ζ y = 0.46 of the pristine borophene to . Our picture provides also a way to decrease the tilt parameter. In this case, the hybridization of the p z orbitals of sites 2, 3 with the R sites (5, 6) reduces as in Tab. I and hence the resultingt and likewiset are scaled down, thereby reducing the tilt to 0.36 in B 6 C 2 -I-[C2&C3]. The above values of the tilt ζ y are calculated based on Eq. (4) and as detailed in SI are in good agreement with the corresponding values directly extracted from DFT bands.
As shown in Fig. S6 of the SI, for doping of two carbon atoms into the structural unit of B 8 , the R-site configurations (c) and (k) have the lower energy than all other configurations, with staggered configuration (k) having slightly lower energy than (c). For a perfect (translationally invariant) doping of two C into the 8P mmn borophene structure, the changes in the tilt are discrete. However, the concentration of C atoms can be continuously varied. Statistical averaging for a concentration x that continuously varies between 0 and 2/8 is expected to generate a continuously varying tilt. The dominant configurations will be the R-site dimers and staggered R-site configurations of panel (c) and (k) of Fig. S6. The lack of perfect backscattering of Dirac electrons is expected to protect the Dirac cone against the Anderson localization.

III. SUMMARY AND OUTLOOK
Based on our ab initio calculations, we have identified the p z orbitals of the I sites as real space sublattice on which the lowenergy degrees of freedom in 8P mmn lattice reside. The effective hoppings between these sites are obtained via renormalization that encodes the virtual hopping via R-sites. This gives a physically clear picture of the formation of the tilt in two-dimensional Dirac cone of 8P mmn borophene. In our picture, tilt arises from a competition between two renormalized hoppingst and t between the second neighbor I-sites. The former involves virtual hopping via two R sites (and 2p x orbitals of the R sites), while the later involves one R site. This builds in a natural difference betweent andt. That is why the pristine borophene has a substantial tilting. Within this picture it is natural to expect that replacing the ridge (inner) atoms by C increase (decreases) the tilt. The rest of the effective hoppings determine the location of the (protected) Dirac node and the anisotropy in the Fermi velocity.
The logic of our work can be applied to other SGs to discover more GQM within a class of 2D materials that afford to provide a "parent" upright Dirac cone [33]: Molecular orbitals in peculiar SGs can mediate longer range hoppings on a backbone lattice of Dirac fermions by promoting it to appropriate graph. The resulting graph in turn deforms the Minkowski spacetime of the Dirac fermions into a metric involving spatio-temporal elements. The ability to tune the tilt parameter is tantamount to controllability of the ensuing solid-state spacetime structure. Dirac materials subject to "gravitational" (i.e. geometric) disorder [34,35] find their salient materialization in the present context when the carbon is randomly substituted for boron. Assuming a doping fraction x of carbon atoms that can be continuously tuned, the most favorable configurations for the C atoms are R-site dimers and staggered R-site ( Fig. S6 (c) and (k), respectively). Random substitutions of this form is expected to lead to continuous variation of tilt as a function of x. The absence of backscattering for Dirac fermions protect it against Anderson localization [36].
The relation between certain graphs and their continuum limit as space geometries is well known [37][38][39]. Our current results suggests a solid-state platform to promote a simple lattice that supports Dirac fermions into a rich graph where vierbeins can be attached to the Dirac fermions [40,41]. This enriches the physics of Dirac materials and promotes them to GQMs where a variety of spacetime geometries can be fabricated that might not even have any analogue in the cosmos. On the 2D materials side, there will be a plenty of room to explore the effects of "geometric" forces or even synthesise non-Abelian gauge fields [42] that are likely to lead to better control in electronic/optical devices, as the effects of spacetime curvature [43,44] can be much stronger GQMs. Furthermore, the geometry of our emergent spacetime in GQMs roots in the Coulomb forces that are much stronger than the gravitational forces. Furthermore, unlike the gravity, the Coulomb forces can be both attractive and repulsive. Therefore the geometric structure of GQMs seems to be much richer than the Einsten's gravity. In three dimensional tilted Dirac/Weyl fermions the resulting geometric theory can be even more interesting: The fermions in such materials are chiral and hence a chiral geometric theory can be relevant to GQM. The developments in GQMs will allow us to emulate aspects of spacetime geometry that is not easily accessible in the cosmos. GQMs have a potential to equip Dirac fermions with types of vielbeins [40] that are beyond the spatial veilbeins of the strain/dislocation paradigm [41].
Supplemental materials to the manuscript: Tunning the tilt of a Dirac cone by atomic manipulations: application to 8Pmmn borophene Since we expect dual readership both from condensed matter and gravity research, in this supplementary information (SI), we provide pedagogical derivations of all the relations and figures referred in the text. To benefit the graduate students, this SI is meant to make the paper self-contained, and understandable on its own trying to minimize consultation with other references.

IV. SYMMETRY ANALYSIS
In this section, we illustrate the application of group theory methods in classification and labeling of the band structures by way of example of 8P mmn borophene. This section is based on chapter 10 of C. Kittel's Quantum Theory of Solids [1]. For a quick and concise introduction to the theory of groups and their representations, we recommend chapter 3 of G. Mahan's Applied Mathematics [2].
Here, the two-dimensional rectangular lattice of pristine borophene B 8 has symmetry group of P mmn. The crystal structure of the pure 8P mmn borophene and the character table for the D 2h (P mmn) point group are presented in Fig. S1(a) and Table II respectively [4,5]. The generators of P mmn space-group are given by two screw is a nonsymmorphic operator meaning the twofold rotation C 2x (C 2y ) around the x(y) axis is followed by a translation a 2 ( b 2 ) as indicated in Fig. S1(a). Note that the inversion center is located at the crossing point of two screw axes (between bond line of atom 1 and 2) as depicted in Fig. S1(a) and the screw is performed around this inversion center. These three independent operations can produce all other symmetry operations of the P mmn space group as follow: Without taking the spin of the electrons into account, D 2h (8P mmn) group has eight 1D irreducible representations given in character Table II. There are two other 2D representations belonging to the double group (i.e. accounting for both orbitals and spin) which is considered in presence of spin-orbit coupling and is not indicated in this table.
First of all, following Kittel [1] we establish that the non-symmorphic elementsC 2y andC 2x that are combination of two-fold rotations around the axis indicated in Fig. S1(a) followed by half-cell translation, guarantee the "sticking together" of the bands along the Brillouine zone (BZ) boundary as indicated in Fig. S1(b). Then we show that band inversion happens by connecting two eigenvalues ± of these operators, as a result of which a Dirac point between Y and Γ is protected by the screw axis operator C 2y , in agreement with our ab initio band structure, as shown in Fig. S1(c).
Suppose that Ψ(x, y) is a solution of the wave equation at k-point along boundary of the BZ (i.e. M X or M Y line). The screw rotationC 2x , mirrorM x , and inversionP operations transform the wave function according toC 2x Ψ(x, y) = Ψ(x + a/2, −y), P Ψ(x, y) = Ψ(−x, −y), andM x Ψ(x, y) = Ψ(−x + a/2, y) respectively which result inM x Ψ(x, y) =C 2xP Ψ(x, y), or the multiplication ruleC 2x =M xP . Other multiplication rules similarly follows. Now to prove that along the M X line defined by k x = π/a, there is a double degeneracy, we show thatP ,M x , and iC 2x satisfy the algebra of Pauli matrices. The definition ofC 2x implies (a is the length of the unit cell in x direction) where, in the the third equality have used the Bloch theorem and in the fourth equality we have used that on M X line we have k x a = π. Subsequently, from the definition ofP andM x , We haveP 2 Ψ(x, y) = Ψ(x, y) andM 2 x Ψ(x, y) = Ψ(x, y). Then by multiplying these operators two by two, we find their anticommunication relations as follows: which indicate that along the M X line defined by k x = π/a, it gives anticommunication relation {P ,C 2x } = 0.
ForP andM x operators, the same procedure yields Naming γ 1 =P , γ 2 =M x , and γ 0 =C 2x , we find thatC 2x ,P , andM x obey the Clifford algebra defined by anticommutation relations: {σ µ , σ ν } = 2Iη µν where I is the identity matrix and η µν = diag(−1, 1, 1) is the Minkowski metric. As is well known the representations of the above algebra are even dimensional, the lowest dimensional representation of which will be two-dimensional where γ µ matrices will be represented by Pauli matrces as e.g. γ 0 = iσ z , γ 1 = σ x , γ 2 = σ y . Therefore along the M X, the bands are at least doubly degenerate. Similar arguments applies to three operatorsM y ,C 2y andP along the M Y direction. Now let us see how the degeneracies can split by moving along the XΓY path. For this we need to label the eigenstates at the Γ point in terms of the eigenvalues ofC 2x(y) . Since the p z states of inner sites in backbone honeycomb sub-lattice are responsible for formation of Dirac cone, the most generic molecular orbitals or Bloch states at the vicinity of Γ point composing the bottom of conduction band |Ψ 1 and |Ψ 4 and those at the top of valence band |Ψ 2 and |Ψ 3 can be constructed from the binding (anti-binding) combinations |p 1 z ± |p 2 z and |p 3 z ± |p 4 z of the 12 and 34 bonds, respectively as, Note that the above ordering that places |Ψ 1 and |Ψ 4 at the bottom of conduction band can be obtained from a simple tight binding (molecular) picture involving atoms 1, 2, 3, 4 in Fig. S1(a). Likewise, it is quite reasonable to place |Ψ 3 (the most "even" combination of atomic orbitals with no minus sign) at the lowest energy, followed by |Ψ 2 at slightly higher energy. Now we apply the action of the point group elements on this molecular orbitals to construct the elementary band representation. As depicted in Fig. S1(a),C 2x replaces the site numbers as 1 ↔ 3 and 2 ↔ 4. Moreover, the |p z orbitals having ∝ z character, change sign under a π rotation around x (as well as π rotation around y). Therefore the effect ofC 2x on the above wave functions becomesC Similarly for the action ofC 2y we have 1 ↔ 2 and 3 ↔ 4, followed by |p z → −|p z (again since p z orbital has ∝ z character). Therefore its effect on the above four states becomes The effect ofP operator on the above wave functions are the same asC 2y . The eigenvalues of operatorsC 2x ,C 2y , andP for bottom of conduction and top of valence band states (+1,+1,+1) and (+1,-1,-1) respectively in agreement with Ref. [3], Now we need all the possible paths where the crystal symmetries connect bands with corresponding eigenvalues of the above three operators in a consistent way. Compatibility relations [1] express that the irrepresentation along the high-symmetry lines are completely determined by the irrepresentation that appear at high-symmetry points (little groups of high-symmetry points/lines). We apply it to the line XΓ and ΓY according to the Eqs. (S11) and (S12). Taking one possible path as an example, as shown in Fig. S1(b) for |Ψ 1 , one of the conduction bands from X point (the eigenfunction ofC 2x with eigenvalue of +1) is connected to the Γ point (the eigenfunction ofP with eigenvalue of +1), and then inverted with the valence band from Y point (the eigenfunction ofC 2y with eigenvalue of +1) that gives (tilted) cone crossing along the ΓY path when we consider the same procedure for |Ψ 2 .
For the formation of cone crossing a "sticking together" (as Kittel puts it) of conduction states and valence states is necessary. Such sticking together is protected by glid/screw elements [1]. In our case it is the screw-symmetryC 2y . Other possibilities are to have the crossing along the ΓX path, or to have a crossing between the blue and red bands of Fig. S1(b). These possibilities for band inversion between states Ψ α with α = 1, 2, 3, 4 is ruled out by the specific ordering given in Eqs. (S9) and (S10).

V. 8PMMN MOTIVATED TIGHT-BINDING MODEL FOR TILTED DIRAC FERMIONS
As pointed out in the main text, effective degrees of freedom are p z orbitals of inner sites the projection of which on the xy plane will be a honeycomb lattice. In this effective model, the buckling of honeycomb lattice sites are ignored. This can be absorbed into effective hybridization with p x orbitals of the ridge sites. Furthermore, the projection of the inner sites of 8P mmn lattice into xy plane will not be a regular honeycomb lattice. Ignoring this will amount to a shift in the location of the Dirac node. This shift can also be absorbed into various effective hoppings. Therefore for modeling purposes, it is enough to start with an "effective regular honeycomb lattice" on which the relevant degrees of freedom (the p z orbitals of the inner sites) reside. In this section we pedagogically derive all the details related to the formation, movement and titling of the Dirac cone based on the effective honeycomb lattice. As is shown in the main text, the picture arising from an effective model can accurately reproduce the DFT tilt parameters. Furthermore, having a tight binding two-band model allows to study various effects such as disorder/interactions/symmetry breaking.

A. Theory without ridge atoms
Without ridge atoms (faint teal sites in Fig. 1(a) and Fig. 2(a) of the main text), life is simple and one basically needs a theory of graphene which will be pedagogically and quickly reviewed below to build upon. We therefore start with Fig. S2 that depicts a regular honeycomb lattice and its associated first Brillouin zone. The three vectors connecting every site from a given sublattice to three nearest neighbors from the opposite sublattice in Fig. S2(a) are , The primitive lattice basis vectors are Using a i . b j = 2πδ ij gives the following basis for the reciprocal space, where a is the distance between two neighbouring atoms. In the hexagonal Brillouin zone it appears the there are six corners K. But in in fact two of them are independent and fall into two categories related to each other by a reciprocal lattice vector constructed from the above vectors. For example the corners of Brillouin zone shown in Fig. S1(b) labeled by K ± a and K ± b are given by, As can be seen adding b 2 to K + a gives K − b and hence they are equivalent. Similarly K − a and K + b are equivalent. So in the calculations, one can pick a convenient one. The effective theories around equivalent K points are the same within a gauge transformation that will be discussed below. The K ± b points indicated in the rhombic Brillouin zone of Fig. S2(b) manifestly emphasizes that there are two independent corners that can not be connected to each other by a reciprocal lattice vectors.
Denoting by a, b the field operators for the annihilation of electrons in p z orbitals of the A and B sublattices of the honeycomb lattice, the parent Hamiltonian is basically the theory of graphene and is given by, Fourier transformation of the above Hamiltonian becomes, that begs for a matrix representation as This is how the spinor structure emerges. The form factor connects the opposite sublattices, and hence is off-diagonal in the above sublattice space and is given by, Plugging in the Cartesian representations of δ 1 , δ 2 , and δ 3 given above and carrying out the summation over the three neighbors, the famous formula of graphene literature can be obtained, The above form factor vanishes at the corners of the hexagonal Brillouin zone. Expanding it around two independent corners, e.g. K ± a gives, The overall factor of −i = e −iπ/2 can be gauge transformed by rotating the overall phase of the wave function in one sublattice. For example gauge transforming the amplitude of wave function on sublattice B as φ B → iφ B , eliminates the −i factor above and we have Upon the above gauge transformation the behavior of the Hamiltonian around the crossing points is given by, Defining 3ta/2 =hv F , this formula gives two representations of the Dirac theory H Dirac (k) = α x k x + α y k y + βm with α x = σ x , α y = ±σ y and m = 0, β = σ z , hence massless Dirac electrons. These two representations are connected to each other by time-reversal operator k → − k followed by T = iσ y K where K is complex conjugation.

B. Role of ridge atoms
The 8P mmn borophene's lattice consists of two different types of B atoms, which we will call inner B (B I ) and ridge B (B R ) atoms in Fig. S3(a). The B I atoms form a hexagonal graphene-like lattice and B R atoms look like one-dimensional chains passing through the effective honeycomb lattice as in Fig. S3(c). As pointed out, in absence of B R atoms, the structure looks like a distorted graphene lattice. This lattice naturally includes the nearest neighbor (black path hoppings in panel (c)) that gives rise to two Dirac cones described above. Deviations of the honeycomb lattice from planar structure allows for mixing with p x orbitals of the ridge sites, while the deviations of the projected honeycomb lattice from perfect C 6 symmetry of graphene accounts for shift in the location of the Dirac cones. Therefore the schematically depicted Dirac cones in panel (d) above are shifted and tilted as described below to form the two Dirac cones in panel (b) of the original 8P mmn lattice.
If there were not ridge sites, in a pure honeycomb lattice, hoppings between 2nd and 3rd neighbors via atomic orbitals would be negligibly small. Therefore there would be no green and red hoppings as in panel (c). However, ridge atoms provide dotted hopping paths of panel (a) that facilitate a hopping between 2nd and 3rd neighbors of the parent honeycomb-like lattice via virtual hopping through the ridge atoms. The connection between the microscopic paths in (a) and effective hoppings in (c) will be discussed in next section. For the purpose of present sub-section, we only need to focus on the effective green and red hoppings in panel (c) above. The corresponding hoppings as depicted in Fig. 2(c) of the main text aret andt for the green path (2nd neighbor). Two different symbols correspond to two different virtual processes depicted in panel (a) above. Similarly for the red paths there are two hoppings t p and t x as depicted in Fig. 2(c) of the main text. As we will see shortly, t p (and t x ) acts like a pseudo-gauge field that shifts the location of the Dirac cone. The superscript x in t x is meant to emphasize a special role played by the p x orbitals of the ridge atoms that will be discussed in the next section. The t p and t x being 3rd neighbor hoppings connect a site from A sublattice to a one in the B sublattice. Therefore the corresponding form factor F xp (k) will contribute off-diagonally to the effective Hamiltonian. Furthermore, the hoppingst andt being 2nd neighbor hoppings connect atoms on the same site, and therefore the corresponding form factorsF (k) andF (k) contribute diagonally to the effective Hamiltonian in the sublattice space. So we end up with where F 0 (k) provides a upright Dirac cone to begin with (as discussed in the previous sub-section).
Using the third nearest-neighbor vectors δ p = − 2 3 ( a 1 + a 2 ) = −2aî, δ 1 where "a" is an "effective" lattice constant for the parent honeycomb lattice, the F xp (k) form factor and its Taylor expansion near the Dirac crossing become, x k y + ia(2t p + t x )k x As can be seen the difference t p −t x between the 3rd neighbor hoppings generated via two different microscopic paths effectively shifts the k y of the opposite valleys by opposite values. Furthermore, the coefficients of the k x and k y above that arise from t x and t p contribute to the re-definition and anisotropy of the Fermi velocity. Remembering to affect the gauge transformation φ B → −iφ B of the previous sub-section, the off-diagonal form factor giving still an upright (but shifted) Dirac cone becomes, This form makes it manifest how the third neighbor hoppings shift the Dirac cones and make their velocities anisotropic. The difference t p − t x is responsible for the mutual movement of the Dirac cones alone the k y direction: The Dirac point at K + a move downward by (t p − t x )/(3at x − 3at 2 ) and Dirac point at K − a move upward by the same value. An illustration of this shift is indicated in Fig. S3(d) by red arrows. So, by increasing difference (t p − t x )/(3at x − 3at 2 ), the Dirac-nodes at K ± a move toward Γ. Equivalently, this means that the Dirac-nodes at K ± b move away from Γ as indicated by red arrows. Formation of the tilt: We are now ready to discuss the green path hoppings in Fig. S3(c) that are second neighbor hoppings on the effective honeycomb structure. The form factorsF (k) andF (k) correspond to the real-space displacements and,δ respectively. TheF Expanding around the Dirac point K ± a , the total AA matrix element becomes, Note that the BB matrix element is the same as above. Remember that the gauge transformation φ B → −iφ B does not change the BB matrix element as a factor (−i) * (−i) gives 1. From the above formula we read: (S15) Eqs. (S13) and (S15) derived for our model system are great guiding principles. Eq. (S13) shows that the 3rd neighbor (red) hoppings are responsible for the location and anisotropy of the upright Dirac cone. Eq. (S15) simply means that the tilt (a long distance property) depends on the difference between the microscopic parameterst −t, namely the difference in the hopping via two different green paths in Fig. S3(c). This model shows how the presence of ridge atoms move and tilt the Dirac cone. In the following section we are going to show how the effective hoppings in panel (c) of Fig. S3(c) arise from microscopic hoppings in panel (a).

VI. EFFECTIVE HOPPING TERMS: RENORMALIZATION VIA MOLECULAR ORBITALS
If one considers the honeycomb-like structure formed by inner sites only, the hopping between 2nd and 3rd neighbors are quite negligible. This is because in a tight-binding approach, the overlap between the atomic orbitals exponentially decays with the distance. But when the ridge sites (denoted by teal color in Fig. S3(a)) are added, they act as virtual sites thorough which a renormalized hopping between the (p z ) atomic orbitals of the inner sites is formed. In this section we provide a simple picture of renormalization based on the molecular orbitals theory that explains why relatively large hoppings between 2nd and 3rd neighbors are obtained. The remarkable agreement of the tilt parameters extracted within such a local quantum chemistry picture indicates its validity. So in this section using molecular orbital theory, we first analytically derive these renormalized hopping parameters of new effective honeycomb lattice in term of hopping parameters of original 8P mmn lattice. This will show that the ridge atoms (B R ) play a crucial role in mediating strong hoppings between the 2nd and 3rd inner sites atoms. Then we discuss the connection between the present molecular orbital treatment with a picture of renormalization procedure by considering Dyson equation for the local quantum clusters of the 8P mmn lattice. The atomic couplings t27 and t37 mix the even parity combination of |2 and |3 with the orbital |7 of the ridge atoms. The resulting molecular orbitals are denoted by red. The odd-parity combination of |2 and |3 is decoupled from the rest, and remains at energy E3 = 0, while a molecular binding orbital constructed from the even combination of |2 and |3 and |7 is stabilized by energy E1 with respect to the initial level 0 of the inner atoms. When the ridge site |7 is eliminated, this lowering becomes an effective hopping between the atomic orbitals |2 and |3 of the third neighbor sites B2 and B3 of the parent honeycomb lattice.
Formation of effective third neighbor hopping: Let us illustrate the basic principle by considering inner sites B 2 and B 3 in Fig. S3(a) and see how the ridge site B7 facilitates a hopping between them. In Fig. S4(a) we have depicted an energy diagram where sites B 2 and B 3 (by symmetry) are at the same atomic energy levels |2 and |3 , respectively. The site B 7 is assumed at state |7 energetically lying 2∆ above the inner sites. The atomic hopping t 27 and t 37 connecting sites B 2 and B 3 to site B 7 that can be derived from Wannier states are indicated. These two (real) hoppings are the same by symmetry. The basic principle can be seen in the limit where t 27 = t 37 2∆. In this limit a virtual second order in ∆ hopping via higher energy atomic orbital |7 stabilizes the energy of the states |2 and |3 by −2t 2 27 /∆. The factor 2 comes from two different ways of performing such a virtual process: |2 → |7 → |3 and |3 → |7 → |2 . The analytical extension of this result to arbitrary 2∆ is straightforward.
In Fig. S4(b) we have denoted the energy level E = 0 and E = 2∆ of the inner and atomic sites before the hybridization t 27 by red and teal colors. The sites B 2 and B 3 are equivalent and hence symmetry adopted basis based on the even and odd representation of a two-element group including identity operation and operator that exchanges B 2 and B 3 is given by, The states |Φ 1 and Φ 2 are even parity and mix with each other by atomic hopping t 27 = t 37 , while the state |Φ 3 being odd parity decouples from the others. In fact state |Φ 1(3) is a bonding (anti-bonding) molecular orbital and is composed of atomic orbitals of inner sites B 2 and B 3 . The matrix elements and Hamiltonian in this basis are As can be seen when t 27 = 0, the bonding/anti-bonding states have zero energy (equal to the energy of atomic orbitals |2 and |3 ) and hence bonding/anti-bonding degrees of freedom are inert in this limit. Upon turning on the atomic hopping to ridge site B 7 , the above structure of the Hamiltonian, still leaves the anti-bonding orbital |Φ 3 decoupled and hence it remains at zero energy, while the bonding (even parity) combination |Φ 1 is mixed with |Φ 2 = |7 giving rise to two split-off states at energies ∆ ± 2t 2 27 + ∆ 2 . The energy of the E 1 = ∆ − 2t 2 27 + ∆ 2 state is always less than energy E = 0 of the original B 2 and B 3 sites. In the low-energy sub-space that -involving only the inner sites B 2 and B 3 -the high-energy site B 7 is eliminated (or integrated out in the renormalization group terminology), this lowering of energy can be interpreted as an effective hopping between B 2 and B 3 , namely, This is so, because in a sub-space composed of only inner orbitals |2 and |3 , placing such a off-diagonal hopping between the two states at energy zero, correctly reproduces the energy lowering of E 1 with respect to E = 0. This effective hopping in the main text has been denoted by t p and is responsible for a pseudo-gauge field that shifts the location of the Dirac node. Note that the in the ∆ t 27 limit, the above espression reduces to the perturbative result −2t 2 27 /∆ cited above, while in the opposite limit ∆ t 27 it becomes − √ 2|t 27 |. Formation of effective second neighbor hopping: As a second example of how longer range hoppings on the parent honeycomb-like lattice of inner atoms are formed, let us consider in Fig. S1(a) the effective hopping between site B 3 in one unit cell and the same B 3 in the lower unit cell (let's call it 3 ). This process will be achieved by two virtual paths: 3 → 5 → 6 → 3 and 3 → 8 → 7 → 3. Let us consider the first virtual path indicated by dotted green path in Fig. S3(a). Denoting by |3 and |3 the atomic orbitals on the two B 3 sites, and by |5 and |6 the atomic orbitals on the intermediate B 5 and B 6 sites, respectively, the four dimensional Hilbert space is broken into even sector and odd sector The hopping matrix elements between the atomic orbitals are shown in Fig. S5 and by symmetry t 53 = t 63 . The hopping between B 5 and B 6 is associated with a π-bonding between the p x orbitals of these two sites. Taking the energy offset 2∆ for the ridge sites with respect to the inner sites into account, the even block of the Hamiltonian and its eigenvalues become, while for the odd sector we have In fact replacement ∆ → ∆ + t 56 /2 maps the odd sector to even sector. This can be interpreted as the fact that the even sector take advantage of the t 56 hopping term between the p x atomic orbitals of the ridge atoms and generate larger stabilization. Therefore for the 3 → 5 → 6 → 3 path, a contribution to the effective second neighbor hopping between the two B 3 inner sites will become, (S17) A similar contribution arises from the path 3 → 8 → 7 → 3 by simply replacing (5, 6) → (8, 7) that gives (S18) Adding these two terms gives the greent in Fig. 2(c) of the main text, (S19) It is rather surprising to note that the Eq. (S17) for the virtual hopping path 3 → 5 → 6 → 3 can be obtained by a simple intuitive picture as follows: As shown in Fig. S5(a), first an antibonding orbital from the ridge sites B 5 and B 6 is formed that has been denoted by φ a . The upward shift of the antibonding orbital with respect to the initial level 2∆ of the ridge atoms is given by the off-diagonal element t 65 in the space of |5 and |6 . Then |φ a = (|5 − |6 )/ √ 2 will have hopping matrix elements t 53 / √ 2 and t 36 / √ 2 to sites 3 and 3, respectively. Then a virtual hopping with these amplitudes via the intermediate state |φ a at energy t 65 + 2∆, precisely gives Eq. (S17). Relation with renormalization: Now we are ready to show that the effective hoppings obtained in this way, admit a nice interpretation in terms of renormalization. As an example, let us elaborate on Eq. (S16) and show that it can be alternatively obtained from a renormalizatioin picture. This will establish that the molecular orbital theory that was explained in the previous part is indeed equivalent to a renormalization procedure.
Imagine the situation we discussed for effective third neighbor hopping energy between B 2 and B 3 . Denoting as before the atomic basis |2 , |3 , and |7 , and breaking the three dimensional Hilbert space into a two dimensional subspace composed of low-energy states |2 , |3 and a one-dimensional high energy subspace described by |7 , the Hamiltonian H can be split as, where the subscripts L (H) stand for low-energy (high-energy) sectors. As described in the textbook of Grosso and Pastori Now we have to evaluate the operator, whose matrix elements are energy-dependent as, The eigenvalues are given by, that immediately becomes, This is exactly the result (S16) of our simple molecular orbital theory. The same logic applies to all other effective hoppings. Therefore the process of virtual hopping via the inner sites that involves molecular orbitals is actually equivalent to formation of renormalized hoppings in the low-energy sector of the theory (residing on the inner sites) as a result of elimination of higherenergy degrees of freedom that reside on ridge sites.

VII. AB INITIO CALCULATION: RELAXATION, ELECTRONIC STRUCTURE, AND STABILITY
The crystal structure of pure 8P mmn borophene is presented in Fig. S1(a). The basic unit cell is rectangular and contains eight B atoms. For DFT calculation we use pseudopotential Quantum Espresso code [9] based on plane wave basis set within the GGA in the Perdew-Burke-Ernzerhof (PBE) parameterization [10]. Simulation of borophene rectangular unit cells is based on the slab model having a 25Å vacuum separating slabs. We also consider monolayer of 8P mmn structure, where some of the B atoms are substituted by C atoms in the form of B 8−x C x (x=0, 1, 2). The obtained structural properties after the ionic relaxations such as lattice parameters, x, y, and z component of the B and C atoms in crystal coordinates are shown in Table III for pure borophene (B 8 ) and C-doped systems B 8−x C x . We have depicted the crystal structure for situation where a single B atom in the ridge (inner) sites is replaced by C atom, denoted by B 7 C 1 -R-[C5] (B 7 C 1 -I-[C2]) in Fig. S6(a) (Fig. S6(b)). An interesting situation where a dimer of B atoms in ridge (inner) sites is replaced by a dimer of C atoms are denoted by B   (1) As we have shown in the main text and the previous sections of SM, the formation of the Dirac nodes is controlled by those atoms residing at the inner sites, while the location and tilting properties are controlled by the atoms residing at the ridge sites. With the aim of finding new materials with tilted Dirac cone, as well as, in order to confirm our analytical model, in this section we investigate whether and to what degree the location and tilt of Dirac-cone would change if boron atoms replaced by another atoms like carbon. We establish that the coarse grained honeycomb lattice of tilted Dirac fermions presented in the main text, is fully supported by ab initio calculations.
For instance, for the case of substitution of single C atom in ridge sites B 7 C 1 -R-[C5], we obtain that lattice parameters and position of atoms are not too different from the pure B 8 (borophene) and, as a consequence, the shape of band structure including location and tilt of Dirac-cone does not change much (see second part of Table III and Fig. S6(e)). On the other hand, the situation is different in the case of substitution of single C atom in inner sites, B 7 C 1 -I-[C2] due to the breaking ofC 2x and C 2y symmetries (it breaks C 2z and inversion symmetries when we consider effective hexagonal lattice). So, single C atoms in hexagonal lattice gap out the tilted Dirac cone bands, as presented in Fig. S6(f). This observation confirms that the atoms at the inner (honeycomb-like) sites are responsible for the formation of the parent Dirac-cone.
In the following, we will consider the situation in which a dimer of B atoms in the ridge (inner) sites is replaced by C atoms B 6 C 2 -R (B 6 C 2 -I). In the case of B 6 C 2 -R-[C5&C6], Table III indicates that atoms in ridge sites get closer to the inner atoms in hexagonal lattice. As a result, the ratio between the effective hopping energy differences and the first neighbor hopping t, namely (t p − t x )/(t − 2t x ) and (t −t)/t (that according to our model in section V of SM, control location and tilting of the Dirac cone, respectively) are expected to be larger than pure borophene. The movement of Dirac cones can be qualitatively seen in panel (g) and (m) of Fig. S6 where a dimer of C atoms are replacing the corresponding B atoms on the ridge sites. A weaker movement can be seen in panel (a) as only one C has been replaced in the ridge site, and hence a weaker shift in the ratio (t p −t x )/(t−2t x ). Fig. S6(m) for B 6 C 2 -R- [C7&C8] indicates that more or less the same story hold when we replace B 7 and B 8 by C atoms. There is another configuration of two C atoms in the ridge sites B 6 C 2 -R-[C6&C7] (see Fig. S6(k)) which is energetically favorable as shown in Fig. S7, but it does not have significant value of tilt. Now let us consider the opposite case, displayed in panel (d) and (h) of Fig. S6 where the C dimers replace B dimers on the inner sites, namely B 6 C 2 -I-[C2&C3] and B 6 C 2 -I- [C1&C4]. In these cases the ridge atoms move away from the xy plane (see Table III), thereby giving rise to opposite effect of panel (c). Therefore in panel (d) the Dirac nodes move away from each other and tilting of the Dirac cone is smaller than pure borophene. This behavior is reinforced in B 4 C 4 -R-[C1-C4] when we replace all inner B atoms by C atoms as the shape of the bands is similar to that of graphene bands (Fig. S6(p)).
So summarize, in terms of the effectiveness of the C substitution in changing the tilt, the most effective substitutions are B 6 C 2 -R-[C5&C6] and B 6 C 2 -I-[C2&C3]. Next priority is the single C-doped. This is because similar to pure borophene they are the C-doped materials in which Dirac cone are dominantly formed by low-energy isolated sets of p z bands. So, in the following we only present the results for pure borophene B 8 , B 6 C 2 -R-[C5&C6], and B 6 C 2 -I-[C2&C3]. The ridge positions in general offer more stable positions for the substitution of C atoms (see Fig. S7). Among the compositions with two C atoms, panel (k) in Fig. S6 has the lowest energy. Then panel (c) has a comparable but slightly higher energy. Among panels (a), (b), panel (a) being ridge site has lower energy.    Due to high flexibility, the range of strain at which 8P mmn borophene remains stable is very much higher than that of other 2D materials such as silicene [14], MoS 2 [15] and black phosphorene [13], and yet, slightly larger than graphene [11] and h-BN [12]. We expect it can also be able to withstand strong strain generated by carbon substitutions. Nevertheless, we need to confirm whether the 8P mmn borophene remains stable by the chemical substitutions with C atoms. Study of phonon dispersion provides a way to investigate the dynamical stability of crystal structures. The calculated phonon dispersions of pure 8P mmn borophene and two systems in which a dimer of C atoms is substituted into borophene (B 6 C 2 -R-[C5&C6] and B 6 C 2 -I-[C2&C3]) are shown in Fig. S8. In this method, the negative value of imaginary phonon frequencies is an indication of dynamical instability. As can be seen in Fig. S8, there are no negative imaginary frequencies, and therefore all three systems studies here are stable with respect to substitution of C atoms for B atoms.

VIII. AB INITIO HOPPING MATRICES ELEMENTS
The parameters of our coarse-grained 8P mmn structure are obtained from atomic scale hopping parameters t ij depicted in panel (a) of Fig. S3. Based on these t ij parameters, as explained in section VI, the effective hopping parameters of our effective model defined on the parent honeycomb structure are calculated. The accurate determination of the effective parameters therefore depends on accurate ab initio calculation of atomic hoppping amplitudes t ij that can be extracted from the corresponding Wannier functions [16][17][18].
The honeycomb lattice tight-binding model presented in our paper to describe the formation of tilt in 8P mmn structure is rather generic, and the hopping parameters can span a wide range of values. Nevertheless for a specific material they are fixed numbers. In contrast to graphene where hopping energy to sites further away than nearest neighbours are much smaller than the first neighbor hopping, our renormalized hopping scenario predict large effective hopping mediated by B R atoms in borophene. Even the graphene fits within our renormalized hopping picture, as in the case of pure graphene, there are no ridge sites, and hence all t ij parameters other than the nearest neighbor hoppings are nearly zero. In contrast, in the 8P mmn structure, they are numerically on the scale of electron volts, giving rise to effective hopping parameters of the same order. So it is sufficient to calculate the atomic t ij parameters. We calculate the hopping energies for all considered systems using Wannier functions. The maximally localized Wannier functions (MLWFs) are constructed with the WANNIER90 library [16,17]. It is worth noting that graphene show well-isolated π bands at the Fermi energy, which induces a simple single band model with p z states. In borophene, our projected bandstructure indicate that the p x states of B R atoms are not completely isolated from the p z states. To verify the adequacy of the above orbitals and hence the validity of calculated Wannier functions, in Fig. S9 for three B 8 , B 6 C 2 -R-[C5&C6], and B 6 C 2 -I-[C2&C3] systems, we compare the DFT-PBE band structure (solid red) with the corresponding Wannier-interpolated bands (blue) obtained with p z Wannier orbitals on the B I site and the p x and p z Wannier orbitals on the B R -site. To avoid complexity, only four bands with p z character are shown. As seen from the band structures, the overall agreement between original and Wannier-interpolated bands is remarkably good. Small deviations appear for states far from the Fermi energy, which given our picture based on the renormalization, are irrelevant for the low-energy physics. Now that we have established in Fig. S9 that the p z orbitals of inner sites and p x and p z orbitals of ridge sites are relevant degrees of freedom, we are ready to calculate the t ij parameters by inclusion of these orbitals. The results of ab initio hopping matrix elements for pure borophene,  Table IV. These values are then used to derive the parameters of the effective model in the right partition of the same table. When a dimer of B R atoms occupying the ridge sites is replaced by a C dimer B 6 C 2 -R-[C5&C6], we obtain comparatively large matrix elements (except t 6,5 ). This is due to the fact that the ridge atoms come closer to the xy plane upon replacement by C atoms as depicted in Fig. 5(a) of main text. So, the ridge atoms mediate stronger hoppings in B 6 C 2 -R-[C5&C6] system, meaning that they can better connect further away inner sites on the effective honeycomb model. As a consequence, the tilt of Dirac-cone is considerably larger than the corresponding values in pristine borophene. The situation is reversed when B I atoms of the inner sites are replaced by a dimer of C atoms B 6 C 2 -I- [C2&C3]. In this case, the ridge atoms move away from the plane of inner atoms, as a results, the hopping parameters are smaller than corresponding ones in pristine borophene. In this situation, the system tends to reduce the tilt of Dirac cone. IV: (left) ab initio atomic hopping matrix elements (in eV ) for borophene and C-doped borophene obtained from Wannier functions. Conventions for labeling of the hopping matrix elements are given in Fig. S3(a). (right) Renormalized parameters of Fig. S3(c). ζy is the tilt parameter calculated from our analytical model. kD/kY quantifies the location of Dirac node extracted from ab initio data. Now that we have a picture of the formation of the Dirac cone, and an analytical model to produce it, we use the model parameters obtained in Table IV to construct a picture of tilted Dirac cone as a function of (k x , k y ) in pure and C-doped borophene. The solution of the tight-binding model leads to the ζ x = 0, ζ y = ±2(t −t)/t, where the further neighbor effective hoppingst andt of the model are calculate in the right partition of Tab. IV. As shown in Table. IV, placing C atoms in the ridge sites generates larger (t −t)/t that ultimately increases the tilt from the ζ y = 0.46 of the pristine borophene to ζ y = 0.59 in B 6 C 2 -R-[C5&C6]. Employing Eqs. (S13) and (S14) the energy dispersion for pristine borophene, B 6 C 2 -I-[C2&C3], and B 6 C 2 -R-[C5&C6] are reconstructed in Fig. S10 using the ab initio effective hopping amplitude reported in Tab. IV. Accuracy of tilt parameters: The ultimate outcome of our model is to understand the mechanism of the formation of the tilt. Hence the value of tilt parameters are important. Let us close this supplementary material by a comparison between the values of the tilt parameter directly extracted from the ab initio data, and the tilt parameter predicted by Eq. (S15) of our model.
We extract the tilt directly from the DFT data by fitting the slopes m R and m L of the ab initio dispersion relations along the ΓY direction as, Using the direct DFT data, we find that the tilting parameter value ζ DFT y are 0.49, 0.47, and 0.66 in pristine borophene, B 6 C 2 -I, and B 6 C 2 -R-[C5&C6] respectively. This is in qualitative agreement with the values obtained from our model, Eq. (S15) and the trend in tilt values upon placing the carbon atoms in ridge/inner sites are correctly reproduced in our model (see Tab. IV). The physical picture emerging from both DFT data and our model is that the substitutions of B R atoms by C dimers where the two C atoms belong to the ridge sites can increase the tilt of Dirac cone. Moreover, the location of tilted Dirac-cone is very sensitive to the position of C atoms on the 8P mmn lattice. This location has been quantified by the parameter k D /k Y to measure the distance from Dirac point to Γ and the results are reported in Tab. IV. If C atoms are replaced in the ridge sites, Dirac-node moves closer to Γ point.
To conclude, the tilted Dirac cone dispersion in 8P mmmn borophene can be well described by our effective tight-binding model. This effective model is not only important for understanding the origin of formation of the tilt in the Dirac cone in borophene, but it also increases considerably the predictive power to find new compounds possessing tilted Dirac cone. Furthermore, out concrete model paves the path for many other studies including the effects of interactions/disorder/symmetry breaking, etc in a physically motivated model of tilted Dirac cone in 8P mmn space group.