Anisotropic scaling for 3D topological models

A proposal to study topological models beyond the standard topological classification and that exhibit breakdown of Lorentz invariance is presented. The focus of the investigation relies on their anisotropic quantum critical behavior. We study anisotropic effects on three-dimensional (3D) topological models, computing their anisotropic correlation length critical exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}ν obtained from numerical calculations of the penetration length of the zero-energy surface states as a function of the distance to the topological quantum critical point. A generalized Weyl semimetal model with broken time-reversal symmetry is introduced and studied using a modified Dirac equation. An approach to characterize topological surface states in topological insulators when applied to Fermi arcs allows to capture the anisotropic critical exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =\nu _{x}/\nu _{z}$$\end{document}θ=νx/νz. We also consider the Hopf insulator model, for which the study of the topological surface states yields unusual values for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}ν and for the dynamic critical exponent z. From an analysis of the energy dispersions, we propose a scaling relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{\bar{\alpha }}z_{\bar{\alpha }}=2q$$\end{document}να¯zα¯=2q and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =\nu _{x}/\nu _{z}=z_{z}/z_{x}$$\end{document}θ=νx/νz=zz/zx for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document}ν and z that only depends on the Hopf insulator Hamiltonian parameters p and q and the axis direction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\alpha }$$\end{document}α¯. An anisotropic quantum hyperscaling relation is also obtained.

Landau's theory of phase transitions 1 provides the standard approach to describe transitions between different states of matter in condensed matter and statistical physics. However, it is not suitable to describe systems that separate phases of matter containing different electronic Bloch states topology [2][3][4][5][6] , also called topological phase transitions (TPTs). In this case there is no symmetry breaking associated with an ordered phase when the system undergoes a TPT 6 . For this reason, many different approaches have been proposed recently to identify the universality classes of topological transitions using scaling ideas 6-18 . In addition, in the late 1990s, it was established the classification table. 1 warning for topological insulators (TIs) and superconductors (SCs) 19,20 that contributed to make this topic a fruitful and exciting area of research. Nevertheless, it is worth to emphasize that there are still many open questions, such as, the fact that some topological phases are beyond this standard classification. For instance, the gapless topological phase and the intrinsic topological (IT) phases 21 . The latter only preserves U(1) charge conservation symmetry and possesses translational symmetry that makes momentum-k a good quantum number 22,23 . IT materials also can exhibit exotic excitations what makes them a very interesting and attractive branch of research [22][23][24][25][26][27][28] .
This research topic is closely related to applications since the discovery of new materials can lead to the development of new devices. Nowadays, there are many technological proposals based on topological materials [29][30][31] . The most prominent applications for TIs rely on their properties of carrying an electric current only at the edges or surfaces in two (2D) and three-dimensional (3D) case respectively, while the bulk remains an insulator 2,4 .
Beyond the standard topological materials, we point out the Dirac semimetals (DSMs) that present band touching at some points in their band structures. These are protected by time-reversal symmetry (TRS) and/ or inversions symmetry (IS) 32 . For this reason, DSMs are recognized to be an intermediate phase between a usual insulator and a TI phase. After a breaking of TRS the Dirac nodes are divided into two nodes, which are called Weyl nodes 33 . In a Weyl semimetal crystal, the chiralities associated with the Weyl nodes (Fermi points) can be understood as topological charges. These leads to monopoles and anti-monopoles of Berry curvature in momentum space, which serve as the topological invariant of this phase 34 . Another property, induced by the projection of the bulk Weyl nodes is the appearance of Fermi arcs in the surface of the Brillouin zone 34,35 . A non-trivial value of the sectorial Chern number C (k z ) guarantees that there are chiral surface states. These form Fermi arcs that connect the projections of two Weyl points with opposite topological charges onto the surface of the Brillouin zone 34 .
Although Dirac semimetals also exhibits Fermi arcs, in the Weyl case these are bulk related while in the Dirac case this is not always true 36 . The Weyl nodes with zero energy observed in Weyl semimetal materials (WSMs) can be described as low-energy fermion excitations using the Dirac equation [37][38][39] .
The band touching in topological semimetals usually exhibits a linear behavior. Nevertheless, Ref. 40 shows that SrSi 2 is a WSM that exhibits a quadratic dispersion. This result inspires us to investigate non-linear band touching, i.e., that exhibits a breakdown of Lorentz invariance 6,7 and consequently is in a different universality class. In this case, the notion of anisotropic scaling becomes fundamental. This is characterized by the ratio θ = ν // /ν ⊥ , where ν // and ν ⊥ denote the correlation length critical exponents along different directions 2,41 .
There is another class of topological materials that is not protected by any discrete symmetry (time reversal, chiral and particle-hole symmetries), represented by IT materials. A special case arises for two band systems, which may realize a Hopf insulator phase. Due to this non-symmetry protected character, a non-trivial Hopf mapping, for m filled and n empty bands m = n = 1 , makes the Grassmannian manifold G m,m+n (C) = G 1,2 (C) topologically equivalent to S 222, 42,43 . From this mapping procedure arises two parameters p and q for which there is a corresponding tight-binding Hamiltonian 22 . If p and q are integers prime to each other (greatest common divisor equal to one) this Hamiltonian describing the Hopf insulator can be explicitly obtained 22,44 .
In this work, we identify different universality classes for anisotropic 3D topological models beyond the standard topological classification, such as, the Weyl semimetal 39 and the Hopf insulator model 22 . We also introduce an generalized Weyl semimetal model to investigate the correlation length exponent ν along any direction. Considering these models, we perform a study of topological surface states to determine the correlation length critical exponent ν , and the dynamic critical exponent z. The former is obtained directly from numerical calculations of the penetration length of the surface states. This length is the characteristic length of the topological phase transition and can be identified as the correlation length. It depends on the distance to the transition and diverges as ξ ∝ |M| −ν where M is the distance to the topological critical point 6,17 .
We use a modified Dirac equation 45 , containing quadratic corrections in momentum, to show that the study of topological surface states is adequate to describe anisotropic features, even for 3D topological models. We show that the only requirement, once there is a non-trivial wave-function solution, is that the surface states only decay into the bulk.
In special, for the Hopf insulator model we investigate the energy dispersion relations and propose a scaling law νᾱzᾱ = 2q , where the critical exponents ν and z depend on only on the parameters p and q and the axis ( ᾱ = x, y, z ) that holds the decaying surface state. We also obtain the anisotropic critical scaling exponent in the form θ = ν x /ν z = z z /z x . We observe that very exotic critical exponents may be possible depending on each pair [p, q]. This is confirmed numerically from the decay of the topological surface state, that for the pair [p, q] = [1,3] yields ν x = 3 and z x = 2.
In general, our results exhibit a connection between the properties of zero-energy surface states and the critical phenomena in the bulk for TPTs beyond the standard classification.

Weyl semimetal model
The Weyl semimetal model with broken time-reversal symmetry T is well-known 34,37,46 . In a cubic lattice, it can be described by the Dirac equation as follows, where t i are the hopping terms, τ i are the Pauli matrices in orbital space and γ is the control parameter of the TPT. Eq. (1) describes two Weyl nodes at positions k 0 = {0, 0, ±k 0 } where k 0 = cos −1 γ with −1 < γ < 1 . When γ = 1 , k 0 = 0 and the two Weyl points merge at a topological phase transition.
Expanding Eq. (1) around the Weyl point k 0 , following Eq. (17) (refer to 'Methods' section for details) we obtain the energy projections in k y and k z directions.
For the projection in k y direction, we take k x = 0 , and at the phase transition point γ = 1 we have This implies a linear crossing behavior of Eq. (2), that is, E(k y ) ∝ t y k y , as we approach the phase transition, since the quartic term k 4 y goes to zero more quickly than the quadratic one k 2 y . The form of dispersion at the transition defines the dynamic exponent z y by E ∝ k z y and we can identify the dynamic exponent z y = 1.
On the other hand, in the k z direction, for k y = 0 , the energy projection becomes Note that, we can rewrite the energy dispersion in Eq.
cos k 0 2 . Accordingly, in a way similar to what happens in the extended SSH model in Ref. 6 , close to the critical point γ = 1 the coefficients are A 1 = 0 while A 2 = 0 . This means that the quadratic term dominates and ensures a quadratic behavior E(k z ) ∝ t z 2 (k z − k 0 ) 2 in the vicinity of the phase transition. This relation also allows to identify the dynamic critical exponent z z = 2.
In Fig. 1 we call attention to this anisotropic band-crossing behavior, linear along the k y (Fig. 1a) and quadratic along k z direction (Fig. 1b). It is well known that this anisotropic way of behaving has impact over the critical exponents along each direction. For this reason, we study the penetration depth of the surface states of the Weyl semimetal model in order to obtain the critical exponents ν and z.
(1) H(k) = t x sin k x τ x + t y sin k y τ y + t z (γ + 2 − cos k x − cos k y − cos k z )τ z where, p = k ( ≡ 1 ), v = 1 , with α i and β are the Dirac matrices, M is the rest mass, B carries the quadratic correction in k and p 2 = p 2 x + p 2 y + p 2 z . In addition, the elements M and B allow to distinguish between the trivial ( MB < 0 ) and non-trivial ( MB > 0 ) topological states 45 .
Next, we put Eq. (1) in the same form of Eq. (4), in terms of γ and expanding directly around the critical From now on, tᾱ > 0 ( ᾱ = x, y, z ) and for this reason B is negative. This means that we have trivial and non-trivial topological phases for M > 0 and M < 0 , respectively. In fact, the non-trivial M < 0 phase region shelters the two Weyl nodes. The point M = 0 , or γ = 1 is the critical point at which the topological transition occurs. The Dirac matrices were replaced by the Pauli matrices in Eq. 5 because the problem has just orbital degrees of freedom and we have two decoupled equations 49,50 .
The investigation of the topological surface states along each direction, requires choosing one direction at a time and considering open-boundary conditions for it. It is important to note that the anisotropy among the momentum directions provides different wave-function behavior as well as different penetration lengths of the surface states in each case.
Surface state conditions. We obtained similar behavior for the wave-functions and penetration depths for open-boundary conditions along x and y directions. This is expected since in Eq. (5) these directions are equivalent. On the other hand, no non-trivial wave-function for the surface state is obtained when z is under open-boundary conditions, see Eq. (34) (refer to 'Methods' section for details). It is reasonable, since this last direction shelters the Weyl points and any projection on a perpendicular plane returns a point.
So, here we consider open-boundary conditions along the x direction (refer to 'Method' section for details). In Eq. (30), we present the proper wave-function and its characteristic decay length ξ ± = −1 ± , see Eq. (26). In order to specify which quantity ξ ± truly represents the penetration length of the surface state, we study the limiting behavior of these lengths in the vicinity of the critical point. From Eq. (26), for γ −→ 1 and k 2 −→ 0 we obtain + → t x B and − → 0 . Since the correlation length must diverge at the phase transition, the length to be identified as the penetration depth or correlation length is ξ − = −1 − −→ ∞. If the surface state only decays in the bulk, the penetration depth ξ − must be a purely real and positive quantity. This means that Eq. (26) must satisfy that returns a narrow interval for k 2 = k 2 y + k 2 z given by Beyond that, the surface states are propagating excitations with energy ε = ±t y k y . . Each blue plane represents the perspective of momentum space with respect to real space, denoted by the gray shape. Along the blue plane, we present the evolution of the surface states in red, pink and yellow colors, as a function of γ . In red color, we have the region where each pair ( k y , k z ) corresponds to a surface state that only decays into the bulk. For this, the condition of Eq. (8) must be satisfied, i.e., − = a is real and positive. In pink color, we have the situation where the surface states decay and oscillate with − = a ± ib , where a > 0 and b > 0 ensure the decay and oscillatory behavior, respectively. The yellow stripes also denote surface states which are purely decaying in the bulk, but distinctively from the region in red these are zero energy states. Accordingly to Fig. 2, as we approach the critical point γ = 1 , the region of surface states decreases until becomes a point. Along this process, we also identify a change of behavior at γ = 0.5 . As γ increases beyond this point, the pink area is no longer observed and just the red region and the yellow stripes, both associated with surface states that only decay into the bulk, remain.
Surface states and critical exponent ν. Notice from Fig. 2 that the only point in the domain of surface states that persists from γ = 0.5 to the critical point γ = 1 is the pair (k y , k z ) = (0, 0) . So, this is the only point that allows to verify the behavior of the penetration depth in the vicinity of criticality.
Accordingly, to perform a study of the topological surface state for this TPT we consider the probability density |ψ| 2 and obtain the site n in the x direction where this has decayed to e −1 of its value at the surface ( n = 1 ) 6 . This allows to obtain the penetration length ξ − as a function of the distance to criticality and consequently to determine the critical exponent ν that controls its divergence. Figure 3, shows the penetration depth ξ − as a function of the distance M to the critical point with openboundary conditions along the x direction. In the plane M − ξ , the behavior of the penetration depth ξ − (red diamonds) as a function of M can be fitted perfectly well by the expression ξ = |M| −ν with the correlation length exponent assuming the value ν = 1 (gray curve). This determines unambiguously that ν = 1 along the x-direction.
We plot the planes ξ − M and |ψ| 2 − x in perspective to show clearly the connection between the divergence of the penetration depth (plane M − ξ ) with the degree of delocalization of the surface state, (plane x − |ψ| 2 ). As we approach the critical point the surface state becomes less delocalized. Then we have identified a diverging length at the topological transition in the Weyl semi-metal, namely the penetration depth of the surface state in the non-trivial topological phase.
solutions where the surface states decay and oscillate with − = a ± ib for a > 0 (pink region), as well as, solutions where the surface states only oscillate (blue region). We also have decaying zero-energy surface states along k y = 0 (yellow stripes). In special, at γ = 0.5 the pink color region vanishes and we observe surface states that only decay. We note, as γ −→ 1 that the radius of the circle containing the surface states decreases until becomes a point at the critical point γ = 1 47,48 .

Generalized Weyl model: emergence of Fermi arcs in all surfaces
In previous sections, we were able to obtain the critical exponents ν x = ν y = 1 from the study of the penetration length of the topological surface states. The dynamic exponents, z x = z y = 1 were, in turn, obtained from an analysis of the dispersion relations at the topological transition, Eq. (1). Since we have just added the term t a sin k z to the τ x in the Weyl Hamiltonian, see Eq. (1), this addition does not break any additional symmetry (refer to 'Methods' section for details). For this reason, Eq. (9) still describes two Weyl nodes but along a different plane, as we will discuss below. In special, we choose t a = √ | 1 − γ | that does not affect the anisotropic character of the energy spectrum and ensures the same critical parameter as before, γ = 1 . This additional term allows us to obtain a non-trivial wave-function in all three directions under open-boundary conditions, including the z-direction. As a consequence, we can obtain Fermi arcs along all surface planes in momentum space, see Fig. 4.
In Fig. 4, we compare the emergence of the Fermi arcs for the Weyl semimetal model, Eq. (1), in Fig. 4a with the Fermi arcs for the generalized Weyl model, Eq. (9), in Fig. 4b. We represent the Weyl points as the blue and red spheres to evidence the positive and negative topological charges, respectively. In Fig. 4a, we can observe the projections of the Weyl points only along the plane k x − k z and k y − k z (four light blue planes), since the Weyl points are along the z direction. On the other hand, in Fig. 4b, we also can see the projections along the plane k x − k y (two light red planes), besides the other planes. This happens because the Weyl points now belong to the plane x − z (for y = 0 inside the cube), and they are not aligned along the z direction, as before. The presence of the term t a sin k z ensures the rotation of the Fermi arcs around the k y direction. We highlight in Fig. 2 the Fermi arcs as yellow stripes, which account for all zero energy, purely decaying surface states.
Next, we perform an analysis of the topological surface states in the generalized Weyl model. This now yields decay exponents ± for open-boundary conditions along any direction (refer to 'Methods' section for details). The main change in ± occurs for open-boundary conditions along x and z directions, where an additional term shows up inside the square root in Eqs. (38) and (39), respectively. In both cases the energy remains ε = ±sgn(B)t y k y and therefore the zero-energy surface states preserve the linear behavior along k y = 0 , as obtained for the Weyl semimetal model. Otherwise, for y open-boundary conditions ± is the same of Eq. (32) and the change occurs at the energy of the surface states now given by ε = ±sgn(B)(t x k x + t a k z ) (refer to 'Methods' section for details). The zero-energy solutions correspond to the tilted Fermi arc along the plane k x − k z in Fig. 4b with linear coefficient −(t x /t a ) . We also verify that these Fermi arcs correspond to the zero-energy surface Physically, the term that distinguishes between Eqs. (1) and (9) restores the penetration depth of the surface state along the z direction with open-boundary conditions. Thereby, the projection of the Weyl points can be found in all planes. The collision between the Weyl points, as the control parameter γ changes, occurs in the plane k x − k z for k y = 0 (into the bulk), that contains simultaneously the red and blue spheres in Fig. 4.
To better understand the trajectory of the Weyl points in plane k x − k z shown in Fig. 4b, we track in Fig. 5 the path followed by each of these points until the critical point at γ = 1 . These points reaarange their trajectories, such that they can collide along the k z direction.
In order to characterize the topological phases and the topological phase transitions, we use the invariant sector Chern numbers 51,52 C (k x ) , C (k y ) and C (k z ) . This invariant is given by C (k µ ) = BZ dk ρ dk F µ with Berry curvature F µ = 1 4π ǫ µ,ν,τ h µ ∂ ρ h ν ∂ h τ /E 3 + (refer to 'Methods' section for details). Here, the indexes ρ and are the orthogonal directions to the µ , ǫ µ,ν,τ is the Levi-Civita antisymmetric tensor and BZ denotes the Brillouin Zone. The topological phases are distinguished by the Chern numbers, C = (C (k x ), C (k y ), C (k z )) . In our case, depending on γ , k x , k y and k z the sector Chern number can be equal to C (k µ ) = 0, −1, 1 , see Fig. 6. Figure 6a,b show three different topological regions, the blue region possesses C (k z ) = 1 (non-trivial), the gray region C (k z ) = −1 (non-trivial) and the white region has C (k z ) = 0 (trivial).
Note that, in Fig. 6 the edges of the region blue in Fig. 6a and the region gray in Fig. 6b correspond to the coordinates of k z and k x , respectively, of the Weyl points for each fixed γ . Accordingly, the projection of these coordinates on the surface delimits the Fermi arc extension, as represented in Fig. 4. Therefore, the topological regions in Fig. 6 are intimately related with the existence of the surface states. For instance, the surface states   Fig. 6a. This solution for the surface states corresponds on the Fermi arcs in Fig. 4 localized at k y = 0 where the plane k x − k z shelters the Weyl points projection (yellow stripe). On the other hand, the surface states with open-boundary conditions only along z direction arise only for k y = 0 and for k x given by the gray region in Fig. 6b. Finally, we also consider open-boundary conditions along y direction to obtain the Fermi arcs numerically. These zero energy solutions are highlighted by the yellow stripe in Fig. 10c (refer to 'Methods' section for details), and depicted in Fig. 4b along the plane k x − k z (light blue plane). In order to study the penetration depth ξ of the zero energy Fermi arcs, we consider the region parameter around the critical point γ = 1 for open-boundary conditions along the anisotropic x and z directions. In Fig. 7, we present the density probability |ψ| 2 obtained applying open-boundary conditions in both directions (middle figure). In cubic perspective, we highlight the density probability for z open-boundary conditions in red color and for x open-boundary conditions in blue color. The decay behavior is drastically different, which is reflected in the penetration depth quantities ξ z and ξ x . In the left figure for z open, we can observe the penetration depth ξ z (red diamonds) as a function of the distance to the critical point M. The exponent ν z = 1 2 gives an excellent fit with ξ z = ξ 0 |M| − 1 2 (black solid line). Otherwise, in the right figure, for x open, the decay of the penetration depth ξ x inside the bulk (blue diamonds) is very well fitted by ξ x = ξ 0 |M| −1 (black dashed line), i.e., with ν x = 1 . So, we are able to identify the anisotropic critical exponent 41 as θ = ν x /ν z = 2 , where x and z are the current anisotropic directions.

Universality class for Hopf insulators: beyond z = 2 dynamic critical exponent
The Hopf topological insulator is a truly 3D topological insulator. It is an example of IT phases, which do not require any other symmetry protection beyond the U(1) charge conservation that is present in all kinds of TIs 22,23 . It is important to note that the Hopf insulator can only be realized in a two band Hamiltonian due to its nontrivial mapping from three-sphere to the two-sphere 22,42 . From the point of view of the Wannier representation, the two-band Hopf insulator can be recognized as a fragile topological system 53 since, simply adding other bands  www.nature.com/scientificreports/ while preserving gap and symmetries trivializes it. On the other hand, in Ref. 54 the authors provide a Hopf multiband generalization that is topological under the addition of the C ′ generalized particle-hole symmetry, such that J −1 HJ = −H * . We also verify the presence of the C ′ for the two-band case, but it is absolutely inherent to the non-trivial mapping that characterizes the Hopf insulator in the first place. Note that, the Hopf insulator remains a no-symmetry protected topological insulator due to the dependence on mapping parameters p and q, as we discuss below. Remarkably, the C ′ symmetry is claimed to protect the surface states in the multi-band generalization, while this is not the case in the two-band model. For instance, consider a general band Hamiltonian with m and n, filled and empty bands, respectively. In this case for m and n different from unit, the space Hamiltonian is mapping into its topologically equivalent Grassmannian manifold and consequently follows the homotopy group Grassmannian classification 43 . For 2 dimensions the homotopy group is π 2 [G m,m+n (C)] = Z and classified by the Chern number. For 3D dimensions and when no additional discrete symmetries are required, the homotopy group is π 3 [G m,m+n (C)] that corresponds only to the identity element. For this reason there is no winding number correspondence in 3D space dimensions 20,55 . However, for the special case of a two band Hamiltonian, i.e., m = n = 1 the topologically equivalent space is S 2 . In this way, for the two band case we have an accidental winding number since the mapping comes from T 3 , a 3D Brillouin Zone torus, to the Grasmannian space parameter S 220, 22 . The non-trivial mapping construction has direct impact on the topological invariant called Hopf index χ 42,56 . In Ref. 22 , the authors show that all Hopf insulators can be realized by the following Here, for each pair [p, q] we have a counterpart tight-binding Hamiltonian in real space which is also able to describe a Hopf insulator model if p and q are prime to each other. We have performed a study of the energy spectra for values of p, and q from [p, q] = [1, 1] to [8,8]. For all cases, we identify the gap-closing points at critical parameter h c = −3, −1, 1, 3 following the high-symmetry points of Eq. (40). The Hopf index is expected to be χ = 0 for |h| > 3 , χ = ±pq for 1 < |h| < 3 and χ = ±2pq for |h| < 1 22 . The non-zero Hopf indexes ensure the presence of gapless surface states in the Hopf insulator phase.
Accordingly, if the dispersion still depends on the critical exponents z and ν , as well as the distance to the critical point M = h − h c we can write for each, ᾱ = x, y, z.
In order to study the gap-closing behavior as M → 0 , we put all k wavevectors at the corresponding highsymmetry points for a specific h c in Eq. (13). This yields E = M νᾱ zᾱ , and we obtain, for all [p,q] values studied, This is a critical exponent relation, independent of h c , as expected. For instance, we get νᾱzᾱ = 2 for [p, q] = [1, 1], [4,1] and νᾱzᾱ = 4 for [p, q] = [1,2], [3,2] , as we can check in Eqs. (41) and (42) for [p, q] = [1,1] and [1,2], respectively. To evaluate each zᾱ , for example z x , we put h = h c , δk y = 0 and δk z = 0 to study the behavior only along x direction, which leads to E ∝ δk z x x . So, in general we expect E ∝ δk zᾱ α for each direction where zᾱ is the dominant term around the critical point. We expand the energy around the critical values of h, for [p, q] = [1, 1] and [1,2] (refer to 'Methods' section for details). In each case we can check the dominant term for each direction and we get z x = 2 for both cases. If we know z x we are able to evaluate ν x from Eq. (14). In this way we obtain ν x = 1 and ν x = 2 for [p, q] = [1, 1] and [1,2], respectively. We perform the same steps for all [p, q] studied here and summarize the results in Fig. 8. The heatmap in Fig. 8 presents the ν x for the corresponding [p, q]. Beyond that, we always find the symmetry ν x = ν y and ν z = 1 . We also observe that the values of ν x follow another critical exponent relation given by, For the present Hopf model, the anisotropic critical exponent 41 can be written as θ = ν x /ν z = 2 for q < p/2 and θ = ν x /ν z = q/p for q ≥ p/2 , since ν z = 1 in both cases. Once again, x and z are the current anisotropic directions.
We also obtain an anisotropic scaling relation for the Hopf insulator, θ = ν x /ν z = z z /z x (refer to 'Methods' section for details), where the x − y plane corresponds to ( ⊥ ) and the z direction to ( ). We derive an anisotropic quantum hyperscaling relation given by, n ↑ = sin k x + it sin k y n ↓ = sin k z + i(cos k x + cos k y + cos k z + h).   www.nature.com/scientificreports/ for details). This relation connects the free energy critical exponent α , to the dimension d, the correlation length critical exponent ν ⊥ , the dynamic exponent z ⊥ and the anisotropic critical exponent θ. Figure 8 gives very unusual values for the correlation length exponents. In order to verify these results, we obtain the penetration depths of the surface states. In Fig. 9a,b we give numerical evidence for ν x = 1 and the more unusual value ν x = 3 , respectively. In order to illustrate the most suitable fit for these cases, we present in Fig. 9 some representative curves for different ν x values using the expression ξ x = ξ 0 M −ν x (label of Fig. 9). For instance, we have ν x = 0.5, 1, 2 and 3 for blue, orange, green and red solid lines. For the case [p, q] = [1, 1] , the best fit in Fig. 9a is for ν x = 1 (orange circles). In a similar way, for the case [p, q] = [1,3] , the best fit in Fig. 9b occurs for ν x = 3 (red circles).
The Hopf surface states are zero-energy modes topologically protected and, as we see, present unusual values for the correlation length critical exponent ν . We recall that this model corresponds to a tight-binding Hamiltonian, as long as p and q are integers prime to each other. The Hopf surface states also penetrate into the bulk, with a characteristic length that diverges at the critical point of the topological transition. In addition, it is expected that a larger Hopf index χ yields a larger number of surface states, but this correspondence remains unclear 42,58 .

Discussion
The main purpose of the present work is to extend and characterize the critical behavior of topological materials that do not belong to the topological table of classification. We have studied two kinds of topological materials, the gapless Weyl semimetal and the no-symmetry protected, intrinsic topological material, the two-band Hopf insulator. The natural spacial anisotropy present in these models allow us to investigate an anisotropic scaling law from energy analysis and the penetration depth of the zero-energy surface states.
The first case we studied of an anisotropic band touching comes from the Weyl semimetal model, where the energy dispersions indicate z y = 1 and z z = 2 for the dynamic critical exponents along y and z directions, respectively. We proposed a modified Dirac equation that better translates this anisotropic aspect to calculate the wave-function of the surface states and its decay exponent. This investigation reveals a trivial, zero wavefunction with open-boundary conditions along the z direction that contains the Weyl points. This precludes the evaluation a critical length exponent ν in this direction.
To overcome this difficulty, we propose a generalized Weyl semimetal model that includes the term t a sin k z τ x to the original model. This addition does not affect the anisotropic behavior, on the contrary, it promotes a rotation of the Weyl points around the k y axis, which enables the presence of Fermi arcs in all planes of the momentum space. The immediate consequence is the possibility to study the penetration depth with open-boundary conditions along all directions. From this, we verify ν x = ν y = 1 and ν z = 1/2 , in agreement with previous results z x = z y = 1 and z z = 2 . We also can determine the anisotropic critical exponent θ = ν x /ν z = 2 . In general, additional terms in the Weyl Hamiltonian H Weyl (k) of Eq. (1) can change the values of the critical exponents, as in the case of different energy dispersion relations at the touching of the bands. However, as long as the anisotropic character of the Hamiltonian is preserved, the scaling law for the anisotropic critical exponent θ = ν x /ν z remains valid. The obtained non-zero sector Chern numbers evidence the non-trivial topological character of these systems. Furthermore, the boundary of region with a non-trivial Chern number in the phase diagram of Fig. 6 is associated with the extension of the Fermi arcs.
In order to go further, we investigate the Hopf insulator model. The non-trivial mapping leads to a tightbinding Hamiltonian for each pair of parameters [p, q], as long as p and q are integers prime to each other. From an analysis of the energy dispersion relations, we were able to determine the critical exponents νᾱ and zᾱ for each ᾱ direction. The values of ν x are presented in Fig. 8 and we discover they follow a pattern summarized in the scaling relation, νᾱzᾱ = 2q . Also we find that the x − y plane is isotropic, such that ν x = ν y = ν ⊥ = 1/2 for q < p/2 and ν x = ν y = ν ⊥ = q/p for q ≥ p/2 , while in the z-direction we always have ν z = ν � = 1 . The anisotropic critical exponent is obtained as θ = ν x /ν z = ν ⊥ /ν � = z z /z x and leads to a quantum hyperscaling relation given by . In this hyperscaling relation, α is the exponent characterizing the singular behavior of the free energy ( f ∝ |δ| 2−α ). Notice, that these results depend only on the Hamiltonian parameters p and q and the anisotropic directions. Indeed, this scaling law holds only for the two-band Hopf insulator, since any additional term may immediately disrupt the very restrict non-trivial mapping that characterizes this model.
Far away from the usual, we obtain exotic values for the critical exponents of the Hopf insulator. For instance, consider p = 1 and q = 3 where the heatmap of Fig. 8 predicts ν x = 3 . In fact, we can check this result obtaining the penetration depth and confirm that ν x = 3 as shown in Fig. 9b. From the proposed scaling relation, we obtain ν x z x = 3 · z x = 6 that renders z x = 2 . Accordingly, using this Hopf scaling law, together with the table of Fig. 8, it is possible to find zᾱ = 2, 4, 6, 8, 10 and even higher values can arise if we extend the table for high values of p and q. In principle, we have not found an upper limit for the Hamiltonian parameters.
We can affirm, based on our results, that even the TPTs beyond the standard classification may support a connection with the theory of quantum critical phenomena. Furthermore, although surface states with non-zeroenergy can exist, only the zero-energy states that are those obeying the bulk-boundary correspondence, allow to make the bridge between the TPT and critical phenomena. The properties of these states are unique and they can yield the set of critical exponents ν and z associated with the spacial and temporal correlations, respectively. To avoid a trivial solution ϕ(x) , we assume that ϕ(x) = χ η φ(x) since τ y χ η = ηχ η with η = ±1 , i.e., χ η is an eigenvector of τ y . Then, Eq. (22) becomes

Methods
The part of the wave-function that depends on x must possess a decay behavior, so φ(x) = e − x where > 0 ensures the decay. Following this, the equation for the characteristic length is given by Then, the solutions for are Since > 0 and assuming that t i ( i = x, y, z ) are always positive, we have η = sgn(B) . Thus, So, we take the general solution for φ(x) that satisfies the Dirichlet boundary conditions, i.e., φ(±∞) = 0 and φ(0) = 0, where ± = ξ ± −1 to recover the usual representation of the penetration depth ξ ± .
The eigenvector χ η for η = sgn(B) is given by The solution for ϕ(y, z) can be found through Eq. (21). Since τ y χ η = ηχ η , we have (16) H(k) = t x k x τ x + t y k y τ y + t z k 2 www.nature.com/scientificreports/ that returns ε = ±sgn(B)t y k y , the energy of surface state that reveals a current along the y direction with velocity v y = ∂ǫ ∂k y = ±sgn(B)t y . So, for the present model along the plane yz we have a current just along the y direction which is consistent with the Fermi arcs results. This also indicates that ϕ(y, z) = ϕ(y) , i.e., it does not depend on z. Thus, the wave-function with open-boundary conditions in x direction assumes the form Surface state wave-function along y direction. In a similar way, for open-boundary conditions along y direction, we obtain the wave-function where, χ η,yo = 1 Surface state wave-function along z direction. Inspired by Eqs. (30) and (31), for open-boundary condition along the z direction, we propose Evaluating the Schrödinger equation for this wave-function at Dirichlet boundary conditions, i.e., H(k)ψ(z, k y , k z ) | z±∞ = 0 where k z = −i∂ z , we have where h x = t x sin k x + t a sin k z , h y = t y sin k y and h z = t z γ + 2 − cos k x − cos k y − cos k z .
For the inversion symmetry, we have H → I H(−k)I −139 where I = τ z that yields So τ z H(−k)τ −1 z = H(k) and the Inversion symmetry is still preserved for the generalized Weyl model. For the time-reversal symmetry, the symmetry operation reads H → T H(−k)T −1 . Since T = IK , where I is the identity matrix and K is the complex conjugate operator, we have Therefore, T H(−k)T −1 � = H(k) and the time-reversal symmetry is still broken. The stability of the two Weyl points requires one of these two symmetries to be broken. If this is not the case, and the material preserves simultaneously I and T we loose this stability and the two Weyl points with opposite topological charges merge leading to a zero total topological charge 39 .
Topological surface state study for the generalized Weyl model. For the generalized Weyl model of Eq. (9) and following the same procedure adopted for open-boundary conditions along x direction, we have the presence of t a k z τ x in Eq. (20). This yields + � = − in Eq. (34). For this reason, we can obtain non-trivial wave-functions for open-boundary conditions along any direction. From this, we can write ± (k y , k z ) in Eq. (38) and ± (k x , k y ) in Eq. (39), that stands for open-boundary conditions along x ad z directions, respectively. For these cases, the surface state energy remains ε = ±sgn(B)t y k y , as for the Weyl semimetal model with openboundary conditions along x. (29) t y k y τ y ψ(x, k y , k z ) = ηt y k y ψ(x, k y , k z )     Accordingly, for open-boundary conditions along y direction, the ± (k x , k z ) preserve the same form obtained for the Weyl model under the same conditions. On the other hand, in this last case the surface state energy assumes the form ε = ±sgn(B)(t x k x + t a k z ) . We also obtain the Fermi arcs numerically for the generalized Weyl model from the diagonalization of the Hamiltonian, Eq. (9), for a lattice with 800 sites, see Fig. 10. We report that the numerical results are in agreement with those obtained via the topological surface state study above. In fact, Fig. 10a,b present the zero-energy ε = ±sgn(B)t y k y = 0 for k y = 0 that describe the Fermi arc as expected. However, in Fig. 10c the zero-energy mode occurs for ε = ±sgn(B)(t x k x + t a k z ) = 0 that yields k z = −(t x /t a )k x , where −(t x /t a ) is the angular coefficient of the linear Fermi arc in the plane k x − k z for t x = 1 , t a = √ |1 − γ | = 1 and γ = 0 . It is important to emphasize, that the extension of the Fermi arcs still depends on the conditions that require that the correspondent must be a real and positive quantity. where, + stands for kᾱ c = 0 and − for kᾱ c = π.