Detecting topological phases in the square-octagon lattice with statistical methods

Electronic systems living on Archimedean lattices such as kagome and square-octagon networks are presently being intensively discussed for the possible realization of topological insulating phases. Coining the most interesting electronic topological states in an unbiased way is however not straightforward due to the large parameter space of possible Hamiltonians. A possible approach to tackle this problem is provided by a recently developed statistical learning method [T. Mertz and R. Valent\'i, Phys. Rev. Research 3, 013132 (2021)], based on the analysis of a large data sets of randomized tight-binding Hamiltonians labeled with a topological index. In this work, we complement this technique by introducing a feature engineering approach which helps identifying polynomial combinations of Hamiltonian parameters that are associated with non-trivial topological states. As a showcase, we employ this method to investigate the possible topological phases that can manifest on the square-octagon lattice, focusing on the case in which the Fermi level of the system lies at a high-order van Hove singularity, in analogy to recent studies of topological phases on the kagome lattice at the van Hove filling.


I. INTRODUCTION
The majority of structural and electronic properties of crystals are ultimately determined by the spatial arrangements of their atoms, which form a lattice structure that repeats periodically in space. In two dimensions, the enumeration of the lattice structures is connected to the problem of finding the possible tessellations of a planar surface, i.e. the different ways to cover an infinite plane by a repeated juxtaposition of certain geometrical shapes (tiles). When taking a single regular polygon as tile, only three possible networks which are homogeneous with respect to vertices, tiles and edges can be formed: the triangular, square and honeycomb lattices [1]. These periodic structures are usually referred to as Platonic lattices and are ubiquitous in condensed matter systems. If the condition of homogeneity is loosened and one is allowed to employ different regular polygons as tiles, the so-called Archimedean lattices can be constructed, which are homogeneous only with respect to vertices [2].
The square-octagon lattice forms one of the eleven possible Archimedean tessellations of the two-dimensional plane [1]. It consists of a repetition of regular square and octagonal tiles, whose vertices define a crystal structure with a four-site unit cell repeated over an underlying square Bravais lattice [see Fig. 1(a)]. The simple tight-binding treatment of the square-octagon nearestneighbor network reveals rather intriguing properties of the electronic band structure: as shown in Fig. 1(b), at 1 /4 and 3 /4 fillings, the energy dispersion shows a partially flat band intersecting two linearly dispersing bands, which form a Dirac cone. The flat dispersion results in the presence of a high-order van Hove singularity [3], namely a power-law divergence of the density of states, * wunderlich@itp.uni-frankfurt.de † ferrari@itp.uni-frankfurt.de ‡ valenti@itp.uni-frankfurt.de which is expected to enhance the effects of electronic correlations and favor the emergence of Fermi surface instabilities [4,5]. From a theoretical perspective, the question of the role of van Hove singularities for the onset of topological phases has been intensively investigated in the recent past in the context of the charge-density wave phase of AV 3 Sb 5 kagome metals (with A=K, Rb, Cs) [6][7][8]. Several works suggested the existence of a topological flux phase among the possible instabilities of the kagome band structure at the van Hove filling [9][10][11][12][13]. It is worth noting that, although a simple tight-binding approach on the kagome lattice yields a band structure with conventional van Hove singularities, recent photoemission experiments detected the presence of a high-order van Hove singularity (close to the Fermi energy) in the band structure of CsV 3 Sb 5 [14]. In this regard, the peculiar electronic dispersion of the square-octagon lattice is an intriguing minimal playground to investigate the possible onset of topological phases when the Fermi level of the system cuts through a high-order van Hove singularity.
It is worth mentioning that there are various proposals of two-dimensional compounds with a square-octagon geometry, such as monolayers of nitrogen group elements [15], metal nitrides and carbides [16], a possible allotrope of monolayer MoS 2 [17], or two-dimensional polymers [18,19]. Most importantly, several synthesis routes to fabricate T-graphene (octagraphene), a tetrasymmetrical carbon allotrope with a square-octagon periodic structure, have been put forward recently [20][21][22][23][24]. The square-octagon lattice is also found as a two-dimensional section of three-dimensional crystals, e.g. in the xzplane of the Hollandite structure, which is characteristic of certain Mn-oxides [25][26][27][28][29][30]. Additionally, the one-fifth depleted square lattice which describes the periodic arrangement of vanadium atoms in the antiferromagnetic CaV 4 O 9 compound [31][32][33] is topologically equivalent to the square-octagon lattice (at first-neighbors) and often referred to as the CaVO lattice. In the past decades, several studies investigated Heisenberg-like models on the  1. (a) Sketch of the square-octagon lattice. The unit cell is delimited by dashed black lines and contains four sites, labeled counterclockwise. Hopping terms are represented by colored lines: first-neighbor bonds in green (t 1 1≤s≤6 ), secondneighbor bonds in blue (t 2 s=1,2 ), third-neighbor bonds in orange (t 3 1≤s≤8 ), fourth-neighbor bonds in red (t 4 1≤s≤12 , dashed lines). For what concerns the latter, only three representative bonds are shown, for the sake of clarity. Arrows illustrate our convention for the complex-valued hopping terms, with j → i denoting the hopping parameter associated with the kinetic process c † i cj. (b) Band structure along the high symmetry path Γ−X−M−Γ of the first-neighbor tight-binding model for the square-octagon lattice, with uniform onsite terms 1≤s≤4 = 1 and uniform hopping parameters t n 1≤s≤6 = −1. The Fermi level has been set at 1 /4 filling. The inset shows the Brillouin zone with its high symmetry points and lines in red.
In this work, we explore the possible topological phases that can manifest on the square-octagon lattice at 1 /4 filling, i.e. when the Fermi level of the system lies precisely at the high-order van Hove singularity. Our study is based on a recently developed statistical learning method [13,51,52], in which a large data set of randomized tight-binding Hamiltonians is generated and subsequently analyzed by statistical tools drawn from machine learning approaches, with the purpose to gain insightful information on possible topological phases. We complement the methodology outlined in Ref. [51] by introducing feature engineering as a tool to identify physical observables that are associated with non-trivial topological phases [52].
The paper is organized as follows: Section II is devoted to the description of the statistical method and discusses the concepts of marginal probability distributions, importance score, and feature engineering, which are employed for the data analysis; in Section III we actualize the statistical method to the specific case of a squareoctagon electronic system, defining the general form of the Hamiltonian; in Section IV we discuss the results of the statistical study, iterating several processes of dimensional reductions of the feature space in order to reach a minimal description of the topological phases; finally, in Section V we summarize our findings.

II. METHOD
Within the framework of the statistical method introduced in Ref. [51], one can explore topological phases on an arbitrary lattice by considering fermionic tightbinding Hamiltonians of the following type The two parts of the Hamiltonian consist of onsite potentials ( i ∈ R) and complex-valued hopping terms (t i,j ∈ C), respectively. The t i,j hopping integrals are taken to be translationally invariant and assumed to vanish when the distance between sites i and j exceeds a certain (arbitrary) threshold. For this reason, we conveniently introduce the notation t n s for the various independent hopping terms, where n runs over all possible Euclidean distances R n , sorted in ascending order, and s is an index denoting the inequivalent bonds at distance n. The electronic filling of the system is fixed by successively filling a certain number of energy bands.
The independent onsite potentials and hopping parameters of the Hamiltonian are dubbed features and collected into the vector x = (x 1 , . . . , x N f ), with N f indicating the number of features. Within the statistical method, a certain choice of the entries of x is referred to as sample and fully determines the electronic properties of the Hamiltonian. At a given filling, each sample can be characterized by the band gap E g , which classifies samples into metals and insulators. If a sample is insulating, we can assign it a topological index, named label, which distinguishes trivial and topological samples. In the following, we choose the first Chern number C as label [53][54][55].
For an unbiased statistical analysis of a particular lattice, we generate a large number of different samples, i.e. different tight-binding Hamiltonian. This can be done by randomly picking tight-binding parameters x i according to a certain probability distribution function (PDF), e.g. a uniform or a Gaussian distribution. The generated data set is then analyzed by calculating the marginal probability distribution functions for each feature x i and label C. Here, ρ C ( x) is the bare PDF of all the insulating samples with Chern number C. Therewith, by inspecting the properties of the marginal PDFs p C (x i ) of the various features x i , we can determine the feature values which are most descriptive for the phase with Chern number C. This allows us to identify, for example, which patterns of hopping parameters is associated to a certain topological phase. We note that features are in general complex-valued. For gaining most insight, one can examine the marginal PDFs for the real part (Re[ The contrast between marginal PDFs for topological phases, i.e. p C =0 , and marginal PDFs for trivial phases, i.e. p 0 , indicates by which features a particular topological phase is characterized. In this regard, a quantitative measure of the importance of a certain feature x i for the topological phase with index C is given by the Bhattacharyya distance [56] between the topological and trivial PDFs, namely This quantity, referred to as importance score in the following, allows to perform a dimensional reduction of the feature space, by omitting the features with lowest values of D B (p C =0 , p 0 ) in the course of the statistical analysis. It is worth mentioning that although other statistical distances between probability distributions can be adopted for the definition of the importance score (e.g., the Hellinger distance [57]), a previous benchmark study on the honeycomb lattice has shown that the Bhattacharyya distance provides a better contrast of the marginal distributions with respect to other metrics [51,52]. Complementary to the use of the importance score, a model can be refined by establishing symmetries between features, as obtained either from physical grounds or from the behavior of the marginal PDFs [13]. An iterative application of the statistical method, involving subsequent data generation, dimensionality reduction and analysis of the marginal PDFs, leads to the definition of effective models for topological phases.
Furthermore, in the present work we pursue a better understanding of the parametrization of topological phases by introducing a feature engineering procedure. We define additional composite features by taking certain combinations, e.g. sums, products or power series, of (some of) the original features and compute their corresponding importance score as the Bhattacharyya distance between the trivial and topological marginal PDFs. Some of these engineered features may carry higher importance score and serve as particularly outstanding descriptors of a particular phase.
Employing the statistical method outlined in this section, complemented by feature engineering, we tackle the study of topological states that can manifest in the square-octagon lattice.

III. LATTICE AND MODEL
The square-octagon lattice, sketched in Fig. 1 (a), is defined by a square Bravais lattice and a unit cell of four lattice sites. Denoting the Bravais lattice vectors by a 1 = (1, 0) and a 2 = (0, 1), the four sites inside the unit cell can be placed at positions ± √ 2 /2 a 1 and ± √ 2 /2 a 2 . To investigate possible topological phases on this lattice we consider a spinless tight-binding Hamiltonian of the form of Eq. (1), with hopping terms being restricted from first to fourth-neighboring sites. Assuming translational invariance, the model contains four onsite potentials, with parameters 1≤s≤4 , and a total of 28 hopping parameters. As shown by the different colored lines in Fig. 1 (a), the 28 hoppings are divided into six first-neighbor terms t 1 1≤s≤6 (green lines), two second-neighbor terms t 2 s=1,2 (blue lines), eight third-neighbor terms t 3 1≤s≤8 (orange lines), and twelve fourth-neighbor terms t 4 1≤s≤12 (dashed red lines; for the sake of clarity, only three symmetryinequivalent links are shown).
The band structure for the model with uniform onsite terms, 1≤s≤4 = 1, and uniform first-neighbors hoppings, Fig. 1 (b). The system is metallic for any filling. The dashed horizontal line indicates the Fermi energy at 1 /4 filling, where the dispersion is characterized by a triply degenerate point at M , formed by the lowest-lying three bands and consisting of a Dirac node and a (partially) flat band. Another band crossing of the same type, formed by the upper three bands, occurs at the Γ point.
By tuning the hopping parameters it is possible to create topological insulators, i.e., open a gap in the energy bands and induce non-zero Chern numbers. In the following, we will infer which parameters have to be manipulated in order to create topological insulators. We focus on the case of 1 /4 filling, for which the Fermi energy intersects the lower triply degenerate band crossing. The presence of partially flat bands gives rise to high-order van Hove singularities in the density of states [5], which implies a strong susceptibility of the system towards symmetry breaking in the presence of electron-electron interactions [3,4]. The Fermi surface of the system coincides with the edges of the Brillouin zone, i.e. it can be seen as the square connecting the M points. Along its vertical (horizontal) edge, i.e. k = (π, k) [k = (k, π)], the Fermi surface displays a mixed sublattice character, with the Bloch waves being evenly localized on sublattice sites 1 and 3 (2 and 4).

IV. STATISTICAL ANALYSIS
For the statistical analysis of the square-octagon lattice we begin by considering the tight-binding model of Eq. (1) with hopping terms up to fourth-neighbor bonds. In order to randomly sample the feature space, we define a set of reference values for each feature, generally denoted as x ref i [51]. The samples are drawn according to a multivariate (two-dimensional) Gaussian distribution in the complex plane, centered in x ref i ∈ C and with covariance matrix Σ = α 2 |x ref i | 2 1 2×2 , where α ∈ R is an arbi- trary hyperparameter. Analogously, for real-valued features, i.e. onsite potentials, a one-dimensional Gaussian is employed. As reference points for the various features, we take ref . Note that the reference points of the hopping terms are scaled by the inverse distance between nth neighbors, i.e. 1/R n . For the width of the Gaussian PDFs, we take α = 0.6. This scheme allows us to consider physical Hamiltonians where extreme values of tight-binding parameters are excluded [51]. For example, within our parametrization, the choice α = 0.6 ensures that the real part of the extracted features does not change sign with respect to the reference point for most samples (≈ 95%). We verified that small changes of α with respect to the above choice do not affect the results significantly. However, in general, extreme values of α shall be avoided. Indeed, if α is too small the sampling is limited to Hamiltonians which are close to the reference point and does not cover a significant amount of the feature space; on the other hand, for a fixed number of samples, choosing a larger value of α leads to noisier marginal PDFs, which may hamper the statistical analysis. We note that the choice of the reference point constitutes the main bias of the present approach. The simplest way to alleviate this bias involves choosing different initial reference points to cover a larger portion of the feature space. The choice can be based on an iterative application of the statistical method: once a set of parameters yielding a certain topological phase is identified, one can perform a new statistical analysis centered around the topological reference point, thus exploring the feature space around it. On the other hand, biasing the results around a certain reference point can be desirable in the case in which the present method is applied to a specific physical system. For example, if one is interested in exploring topological phases for a certain target material, the reference point can be chosen to be an ab-initio determined tight-binding Hamiltonian [52].
After creating a data set of n S = 2 · 10 7 samples on the square-octagon lattice, we find 3.3% insulators out of all samples, 17.6% of which are topological. Nearly all topological insulators (99.6%) have Chern index C = ±1.
As we sample in a large parameter space, the number of topologically non-trivial samples is small. Hence, we proceed attempting a dimensional reduction in order to infer more information on the topological phases.

A. Dimensional reduction
The parameter space can be reduced by examining the importance scores D B (p 1 (x i ), p 0 (x i )) and D B (p −1 (x i ), p 0 (x i )) for the C = ±1 phases, which constitute the majority of topological samples. The importance score of each feature is the same for C = 1 and C = −1, because the underlying marginal PDFs show specular behavior with respect to the Im(x i ) = 0 axis in complex plane for opposite Chern numbers. As shown in Fig. 2, we observe zero importance for onsite terms. Hence, the parameters i do not play any role in distinguishing trivial and topological phases, and can thus be excluded from the statistical analysis. On the other hand, the importance score of all hopping parameters is finite.
Similar values of the importance scores, which vary up to statistical noise due to the finite sample count of topological insulators, indicate the presence of sub-groups of hoppings, as expected from the inherent symmetry of the square-octagon lattice. Indeed, we can distinguish two classes of first-neighbor hoppings, according to their importance score: (i) bonds within square plaquettes {t 1 1 , t 1 2 , t 1 3 , t 1 4 } and (ii) bonds connecting square plaquettes {t 1 5 , t 1 6 }. Also fourth neighbor hoppings can be grouped in three classes: (i) bonds crossing the square plaquettes {t 4 1 , t 4 2 , t 4 3 , t 4 4 } [vertical red dashed line in Fig. 1 (a)], (ii) bonds crossing the octagonal plaquettes and connecting sites belonging to the same sublattice {t 4 5 , t 4 6 , t 4 7 , t 4 8 } [horizontal red dashed line in Fig. 1 (a)] and (iii) bonds crossing the octagonal plaquettes and connecting sites belonging to different sublattices {t 4 9 , t 4 10 , t 4 11 , t 4 12 } [diagonal red dashed line in Fig. 1 (a)]. For what concerns second-neighbor hoppings (t 2 s=1,2 ) and third-neighbor hoppings (t 3 1≤s≤8 ) no distinction into sub-groups can be made based on the importance score.
Among all hoppings, the lowest importance is shown by the fourth-neighbor hoppings which connect sites belonging to the same sublattice, i.e. classes (i) and (ii). Therefore, based on this observation, we omit these parameters (and the onsite potentials) in the next iteration of our analysis. A new sampling procedure with the reduced model yields a remarkably larger number of insulators (64% of all samples) and shows a rather low importance for the second-neighbor hopping terms, which is approximately three times lower than the importance of third-and fourth-neighbor terms (not shown). Hence, based on this observation, we also exclude the secondneighbor hoppings from our model, in order to scale down the size of the feature space. This will enhance the contrast between the marginals PDFs and thus simplify the subsequent analysis.
We are thus left with a model including only first-, third-, and fourth-neighbor hoppings of class (iii) (i.e., those connecting sites belonging to different sublattices). Note that, for simplicity, we will refer to the latter as "fourth-neighbor hoppings" in the remainder of the paper. The new data set contains 64% insulators, 18.2% of which possess a non-trivial Chern index C = 1 or C = −1. The fraction of insulators with higher Chern number is negligibly small. Compared to the previous iterations, we observe a higher portion of topological insulators due to the reduced parameter space (11.6% out of all samples, against 0.57% for the full model including onsite terms and all hoppings up to fourth-neighbors).
We can gather information on topological phases from this model by considering the marginal probability distributions. Based on their appearance, the PDFs of firstneighbor hoppings can be grouped in two subsets, one formed by the hoppings inside the square plaquettes {t 1 1 , t 1 2 , t 1 3 , t 1 4 }, and the other containing hoppings that connect adjacent plaquettes {t 1 5 , t 1 6 }, see Fig. 1 (a). The marginal distributions for third-and fourth-neighbors, respectively, show the same behavior among each type. One exemplary set of the PDFs of imaginary parts p C (Im[t n s ]), which indicate the "directions" of the complex hoppings, and PDFs of the moduli p C (|t n s |), which describe the overall hopping strengths, is shown in Fig. 3 for each group of hoppings. From these PDFs we can infer the most descriptive features characterizing trivial and topological phases, as discussed in the following.

Trivial C = 0 phase
In the trivial phase, the marginal PDFs for the imaginary parts of all hoppings shown in Fig. 3, i.e., p 0 (Im[t 1 1 ]), p 0 (Im[t 1 5 ]), p 0 (Im[t 3 1 ]) and p 0 (Im[t 4 5 ]) show a perfect symmetric behavior around zero. We can thus infer that no specific hopping direction is preferred. The PDFs of the moduli p 0 (|t 1 1 |), p 0 (|t 1 5 |), p 0 (|t 3 1 |) and p 0 (|t 4 5 |) show similar shapes, with a non-zero mean indicating finite bond strengths. Hence, the trivial insulating phase can be realized by finite first-, third-and/or fourth-neighbor hoppings, with no specific hopping direction (e.g., by real hoppings). This configuration is schematically illustrated in the left panel of Fig. 4, where we color the relevant bonds within one unit cell.

Topological C = ±1 phases
For first-neighbor hoppings within the square plaquettes, exemplified by the term t 1 1 in Fig. 3, we observe that the marginal PDFs for non-zero Chern numbers, i.e., p ±1 (Im[t 1 1 ]) and p ±1 (|t 1 1 |), look rather distinct from the PDFs of the trivial phase. For the C = 1 phase, Im[t 1 1 ] tends to be larger than zero which corresponds to a counter-clockwise winding of the hoppings around the square plaquettes. p −1 (Im[t 1 1 ]) is the conjugate of p 1 (Im[t 1 1 ]), hence the winding is clockwise. The modulus |t 1 1 |, i.e., the overall hopping strength, shows larger values for topological phases than for the trivial phase. For what concerns the remaining first-neighbor hoppings, represented by the term t 1 5 in Fig. 3, we observe that p ±1 (Im[t 1 5 ]) is symmetric around zero, i.e., no particular hopping direction is indicated and, thus, these hoppings do not play a role in differentiating between C = ±1 and C = 0 phases. At variance with the case of |t 1 1 |, the marginal PDFs of |t 1 5 | have similar means for C = ±1 and C = 0 phases.
As shown in Fig. 3 by the representative term t 3 1 , also the marginal PDFs for the imaginary part of thirdneighbor hoppings behave differently for topological and trivial phases: Im[t 3 1 ] tends to be larger than zero for C = 1, while for C = −1 it shows a tendency to be smaller than zero. This implies that non-zero thirdneighbor bonds with complex hoppings winding clockwise (anti-clockwise) in the octagonal plaquettes can support the non-trivial C = 1 (C = −1) phase. On the contrary, the marginal PDFs of the moduli p C (|t 3 1 |) look identical for C = ±1 and C = 0 and, thus, they provide no information about the topological properties. Finally, we observe that the marginals of fourth-neighbor hoppings, represented by t 4 5 , show qualitatively the same behavior as the marginals for the third-neighbor bonds. Hence, also non-zero fourth-neighbor bonds with hopping directions winding clockwise (anti-clockwise) can support the non-trivial C = 1 (C = −1) phase.
In summary, a topologically insulating phase with C = 1 can be induced by anti-clockwise first neighbor hoppings on the square plaquettes, which are relatively stronger than the bonds connecting square plaquettes, together with clockwise third and fourth neighbor hoppings. Topological insulators with C = −1 can be created by reversing the hopping directions of the C = 1 phase. These results are schematically summarized by the sketches in the middle and right panel of Fig. 4. Here, the thickness of the bonds reflects the relative hopping strengths and the arrows illustrate the hoppings directions, i.e., the sign of their imaginary parts.

B. Towards a first-neighbor model and feature engineering
To gain a deeper understanding of the phases which can manifest in the square-octagon lattice, we continue with a reduction of parameters based on the importance scores. Within the tight-binding model with first-, thirdand fourth-neighbor hoppings discussed in the previous section, the importance score of first-neighbor terms turns out to be up to eight times larger than the importance score of third-and fourth-neighbor terms. Based on these observations we exclude third and fourth neighbor terms as the next step of our analysis.
Creating a data set for the Hamiltonian with only firstneighbor hoppings yields 95.5% insulating samples. The fraction of topological insulators corresponds to 13.7% of all samples, analogously to what has been observed in the calculation with first-, third-and fourth-neighbor hoppings. This further indicates the higher importance of first-neighbor hoppings for the topologically non-trivial phases with C = ±1. With the first-neighbor model we arrive at the minimal possible description for topological phases on the square-octagon lattice. As done previously, we can group first-neighbor hoppings into subsets based on the behavior of marginal PDFs (not shown): (i) hoppings that form square plaquettes, {t 1 1 , t 1 2 , t 1 3 , t 1 4 }, and (ii) hoppings connecting different squares {t 1 5 , t 1 6 }. In order to try gaining additional information on the topological phases, we apply feature engineering, namely we define new composite features by taking all possible products involving (distinct) first-neighbor hoppings, i.e. pair-wise products of the form t 1 s t 1 s , triple products of the form t 1 s t 1 s t 1 s , and so on, up to the product of all six firstneighbor hoppings. We then calculate the importance score for the newly engineered features and identify the ones which play a major role in characterizing the topological phases. As shown in Fig. 5, the product of all hoppings on the square plaquettes, namely t 1 1 t 1 2 t 1 3 t 1 4 , turns out to possess a remarkably large importance score, D B = 0.45 (c.f. D B ≈ 0.05 for t 1 1≤s≤6 ). For this particular engineered feature, the marginal PDFs in the complex plane, shown in Fig. 6, provide crucial insight. In the C = 0 phase, the PDF of t 1 t 2 t 3 t 4 is symmetric with respect to the real axis, as shown in Fig. 6 (a). On the other hand, the marginals for the topological phases [ Fig. 6 (b) and (c)] are completely localized in the upper and lower half of the complex plane for C = −1 and C = 1, respectively. Hence, the importance score for distinguishing the two topological phases takes its maximal value, i.e., D B (p 1 (t 1 1 t 1 2 t 1 3 t 1 4 ), p −1 (t 1 1 t 1 2 t 1 3 t 1 4 )) = ∞. This implies that the distinct topological phases are unambiguously distinguished by this engineered feature. Physically, the topological phases are distinguished by the phase picked up after one loop in the square plaquette which is given by ϕ[t 1 1 t 1 2 t 1 3 t 1 4 ] = ϕ[t 1 1 ] + ϕ[t 1 2 ] + ϕ[t 1 3 ] + ϕ[t 1 4 ]. Eventually, the engineered feature t 1 1 t 1 2 t 1 3 t 1 4 may serve as the unique descriptor of the topological phases.

V. SUMMARY
The statistical method introduced in Ref. [51] constitutes an effective procedure to identify possible topological phases that can be realized by a tight-binding Hamiltonian on a given lattice structure. We employed this technique to scrutizine the topological phases appearing at the high-order van Hove filling on the square-octagon lattice, which forms one of the eleven Archimedean tessellations of the two-dimensional Euclidean plane and is realized in a number of different materials. Starting from a generic tight-binding model with hoppings up to fourth nearest neighbors, we constructed a dataset of randomized Hamiltonians labelled by their Chern number as topological index. We then performed a statistical analysis of the marginal probability distributions for the various parameters of the system and, by means of dimensional reduction, we reached an effective model describing topological phases with Chern number C = ±1 on the square-octagon lattice. Most importantly, we introduced a feature engineering procedure that allows us to gain deeper insight into the nature of the topological phases by identifying polynomial combinations of tight-binding parameters which are associated to non-trivial topology, e.g. Peierls-like fluxes. Going beyond the methodological improvements presented in this work and the results for the square-octagon lattice, the present statistical method can be regarded as a potential tool to perform a materialspecific search of topological phases, by exploring the phase space around a tight-binding Hamiltonian obtained from first principles (e.g. by density-functional theory and Wannierization). The search of topological phases and its characterization by means of engineered features could serve as a guide to experimental manipulation of the target material to tune its properties towards desired topological phases, e.g. by means of applied pressure or strain.