{'fraction_non_alphanumeric': 0.030422944273007804, 'fraction_numerical': 0.01383191141768016, 'mean_word_length': 4.5024970035956855, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evolving Heterotic Gauge Backgrounds: Genetic Algorithms versus Reinforcement Learning

The immensity of the string landscape and the difficulty of identifying solutions that match the observed features of particle physics have raised serious questions about the predictive power of string theory. Modern methods of optimisation and search can, however, significantly improve the prospects of constructing the standard model in string theory. In this paper we scrutinise a corner of the heterotic string landscape consisting of compactifications on Calabi-Yau three-folds with monad bundles and show that genetic algorithms can be successfully used to generate anomaly-free supersymmetric SO(10) GUTs with three families of fermions that have the right ingredients to accommodate the standard model. We compare this method with reinforcement learning and find that the two methods have similar efficacy but somewhat complementary characteristics.


Introduction
Obtaining string theory solutions that match the properties of particle physics in our world has been one of the main driving forces in string phenomenology, fuelling a significant effort in developing geometric engineering tools for constructing physically viable string models. This effort has met with only partial success due to the sheer mathematical difficulty of the problem and, to date, none of the existing string theory models allows for a proper description of particle physics. However, the model building experience of the past few decades has taught us much about the magnitude of the problem and about how (and also about how not) to approach it.
The first lesson we have learnt is that the size of the string landscape is much larger than previously thought. The famous first estimate of O(10 500 ) consistent type IIB flux compactifications [1,2] seems rather conservative in comparison with the latest estimates. For instance, in Ref. [3] it was shown that a single elliptically fibered four-fold gives rise to O(10 272,000 ) F-theory flux compactifications. The second lesson is that the number of compactifications that match the symmetry group and the particle spectrum of the Standard Model is very large, despite representing only a tiny fraction of all consistent compactifications to four dimensions. In Ref. [4] it was argued that there are at least 10 23 and very likely up to 10 723 heterotic MSSMs, while the authors of Ref. [5] argued for the existence of a quadrillion standard models from F-theory.
The third lesson is that these numbers are so large that traditional scanning methods cannot be used for systematic exploration. One can, of course, focus on small, accessible corners of the string landscape and this approach has been successful to some extent. For instance, in Refs. [4,6], 10 7 pairs of Calabi-Yau three-folds and holomorphic bundles leading to SU (5) heterotic string models that can accommodate the correct MSSM spectrum were explicitly constructed. However, going beyond just the spectrum is both computationally challenging and may well lead to an empty set of models, depending on the required degree of similarity with the standard model. What is needed instead is a tool-set of non-systematic search methods, that can quickly detect phenomenologically rich patches of the string landscape and which can implement 'on the fly' a number of checks that go beyond the usual spectrum considerations.
Such methods can be developed using modern techniques of optimisation and search, and have already been implemented in a number of string contexts. Concretely, genetic algorithms (GAs) were introduced as a search method for finding string vacua with viable phenomenological properties in Ref. [7], with initial applications in the free fermionic formulation of heterotic string theory. More recently, reinforcement learning (RL) has been used in Ref. [8] to generate type IIA intersecting brane configurations that lead to standard-like models and in Refs. [9,10] to identify SU (5) and SO(10) string GUT models, respectively, that can lead to three generation MSSMs. The landscape of type IIB flux vacua was explored in Ref. [11] using GAs and Markov chain Monte Carlo methods, while in Ref. [12] the same methods were used, as well as RL.
In the present paper we use a GA to obtain phenomenologically viable heterotic string compactifications on two Calabi-Yau three-folds with monad bundles. The same class of models has recently been explored in Ref. [10] using RL, and indeed one of the purposes of the present work is to compare the two methods. We search and obtain several hundred SO(10) models that can accommodate three standard model families of quarks and leptons. We also check during the course of the genetic evolution that the models are anomaly free and that they pass a necessary condition for bundle stability. The latter check would normally be performed at the end of the search, but recent advances in establishing simple, analytic formulae for line bundle cohomology [13][14][15][16][17][18][19][20] have made possible the inclusion of additional checks during the search.
The rest of the paper is organised as follows. In Section 2 we give a brief introduction to GAs while in Section 3 we review the construction of heterotic models on Calabi-Yau three-folds with monad bundles. In Section 4 we discuss the details of our GA implementation and in Section 5 we present the results. 2 The GA technique for seeking perfect models Let us begin by giving some background and motivation for the GA. The approach of GAs (and evolutionary computation more generally) is to use concepts that are generally understood to underly evolution in order to find an optimum solution to a specified problem. Thus all GAs contain an evolving population of individual solution attempts, a "fitness" function for each individual, and a "selection and breeding process" based on their fitnesses, which applies evolutionary pressure with the aim of driving the population towards a solution. Depending on the problem at hand, success may be declared when the population contains sufficiently many "good" examples, or when one or several individuals in the population satisfies every one of a discrete set of desiderata. The present study will be of the latter kind. More specifically, we will be searching for models with properties that make them candidates for a standard model from string theory. For easy of terminology we will refer to such models as "perfect models" in the following.
The idea of evolutionary computing was developed in the 1950s, but it was formalised by Holland [21] who moulded the idea of "building blocks" into a theory of "schema". Later reviews can be found in Refs. [22][23][24][25][26][27][28][29][30]. There has been growing interest in the technique within the astrophysics and particle theory community, for example in Refs. [7,11,, partly due to its often surprising efficacy.
We will return to discuss some of the ideas that have been proposed to explain the efficacy of GAs in a moment, but let us first describe the technique itself in more detail as it pertains to the problem under discussion. The fundamental starting point of any GA is to encode the problem into a string of data, referred to as the genotype (each entry of which is referred to as an allele).
To be specific, in the present case the genotype is constructed as follows. Fixing a Calabi-Yau three-fold X, we seek monad bundles V constructed from two sums of holomorphic line bundles B and C, with ranks r B and r C , respectively, via the short exact sequence 0 → V → B → C → 0. Each line bundle in B and C is specified by its first Chern class, which, relative to a basis of H 2 (X), is described by an integer list with length equal to the Picard number h = h 1,1 (X) of X. This means that V is specified by h(r B + r C ) integers. In fact, as explained below, we will be looking for bundles satisfying c 1 (V ) = 0, a constraint that cuts down the number of integer parameters to h(r B + r C − 1). In order to translate this into a genotype we will allow the line bundle integers to assume values only within a finite range, which can then be implemented as a binary encoding, with each integer being mapped to a binary string of fixed length n B for B and n C for C.
Thus the collection of all line bundle integers in B and C or, rather, the corresponding bit-string obtained by joining their binary encodings defines a particular model. The starting point of the GA is then a population of N pop randomised individuals, each being represented as a string of digit alleles of length = h(r B n B + (r C − 1)n C ) .
(2.1) For instance, for the first of our working examples we will have h = 2, n B = n C = 3, r B = 6 and r C = 2, giving = 42.
In the case of such a binary encoding there are general arguments [22,25,28] that suggest that the optimum size of the population N pop depends only logarithmically on the genotype length , with a formal lower bound determined by requiring that both possible values for each allele should appear in the initially randomised population with some representative probability P * . (For example P * = 0.999, would imply that there is a 1 in 10 3 probability that any individual allele contains only a 0 or a 1 in the entire population.) Typically this gives an estimate of N pop 20 that is a significant underestimate of the optimal size of the population, but it nevertheless gives an indication of the required population size. For the current study, the GA was found to perform well with N pop 100. The GA then consists of repeated application of the following basic elements on the population: Selection: Individuals are first selected from the population to make N pop "breeding pairs" which will create a new generation to replace the original one. In order to do this the first step is to calculate the properties of each individual (its so-called phenotype). This particular step of extracting the phenotypes from all the genotypes is the true hallmark of a GA, as the relation between the genotype and the phenotype it generates, and also the correlations between different properties within the phenotype, may be very deeply hidden. This step, which is where the physics comes into play, is also often the most case-dependent and computationally-intensive part of the whole procedure. Encouragingly, however, it is the aspect that is most amenable to parallelisation.
Based on the phenotype, each individual is then assigned a fitness value, which is a function of the individual's phenotype designed such that it increases as the model approaches a desired optimal solution. There is a wide variety of methods for assigning the fitness to an individual, and the process may take many different forms, such as rank-weighting, tournament selection, and so on. Designing a good fitness function is a crucial factor in the success (or otherwise) of the GA, because its convergence will depend strongly on the correlation between the fitness and the distance from the hypothetical optimal solution in genotype space (that is, the Hamming distance). This fitnessdistance-correlation was formalised in a number of works by Jones, Forrest et al [24,[52][53][54][55].
We can appreciate this aspect of the GA for the simple example of maximising a smooth function. Suppose we choose a fitness functions which is simply proportional to the function value. Then, as the population approaches a relative flat maximum, there may be little difference in fitness between individuals which are nevertheless widely separated in coordinate space (which in the case of a continuous function would be directly given by the genotype). In other words there would be little fitness-distance-correlation in the neighbourhood of a solution and convergence would cease. A simple way to avoid this kind of convergence problem is to rank the individuals in the population by their fitness value, and then to simply equate their fitness with their ranking. This is the aforementioned rank-weighting selection procedure, and is the approach that we will adopt in this study. It has the additional advantage that we do not need to store all fitness values separately, but can keep track of the ranking by re-ordering all the N pop genotypes within a single N pop × array. Selection to participate in breeding pairs is then carried out with a linearly distributed likelihood (known generally as roulette wheel selection) that simply depends on an individual's ranking. In the present case an individual at rank k in the population is selected for breeding with a probability that depends linearly on its ranking as where α determines how much more frequently the fittest individual is chosen over the least fit. In other words, the probabilities P Npop and P 1 for selecting the fittest and least fittest individual are related by Typically, an optimal configuration is found when the fittest individual is selected for breeding a few times more than the least fit, so that a good choice would be α = 2 or 3. But note that it is important that less fit individuals are also able to breed.
Breeding/Cross-over: At this stage we have N pop breeding pairs, with the median individual having been selected for breeding twice (that is, N pop P Npop/2 = 1 + O(1/N pop )). A new population of individuals is then formed by splicing together the genotypes of the two individuals in each breeding pair. Again there are many different ways to do this. The canonical choice is N -point cross-over, namely one cuts the genotypes of a breeding pair at the same N random positions along their lengths and swops the cut sections. For robustness, the publicly available codes often use sophisticated combinations (for example [26,31,56,57] uses both one-and two-point cross-over in roughly equal proportions) to reduce end-point biasing. In the present study this effect was not found to be significant, and a simple one-point cross over was found to work sufficiently well.
Mutation: With only the two previous elements, one would already observe clustering of the population around good solutions over generations. However, the power of GAs comes from the third element which is mutation. Once a new generation is formed, a small fraction (usually around a percent) of the alleles are flipped at random. This prevents stagnation in the population, where the entire population clusters around a local maximum in the fitness, when better solutions are available globally. For example, it removes the possibility alluded to below Eq. (2.2) that a particular allele could end up being blocked with either all 0's or all 1's running throughout the entire population.
It is important to understand that mutation is not just an improvement to the convergence, but is absolutely integral to the entire process. Depending on the problem and the structure of the fitness landscape, it causes a dramatic increase in the overall rate of convergence. This can be seen in practice by optimising the mutation rate as demonstrated in Ref. [7]. It is worth mentioning one possible innovation that could be considered, but which we did not use in this study, namely so-called creep mutation. This is effectively a small mutation in the phenotype which overcomes socalled Hamming walls, which occur when the population is close to an optimum solution in terms of phenotype, but far away in terms of Hamming distance. An example would be the binary number 1000 and 0111 which are numerically close but have four different digits. This can again be thought of as a way to overcome breakdowns in the fitness-distance-correlation, and could be thought of as assisted mutation.
Elitism: Finally to create the new generation in a way that guarantees a monotonically increasing maximum fitness, we invoke elitism. This means the fittest individual from the previous generation is copied into the new one and replaces the least fit new individual.
The process is then repeated over many generations, and is terminated and reset when each "perfect model" is found. Note that there are many other practical elements that could be considered, such as fitness "crowding penalties", and "niching", that we do not discuss (or require). They are covered in the literature (see Refs. [25,27]). The overall procedure that we use in this study is summarised in the following pseudo-code:

Discussion of efficacy
Before moving on to the details of heterotic string models with monad bundles, it is worth pausing to consider why GAs are thought to work. With the benefit of several intervening decades of development in heuristic search methods, we can see that GAs contain an interesting mix of elements that can be found in other methods. For example the fitness function is arguably just a costfunction by another name. In addition, the general idea of fitness-distance-correlation in Jones, Forrest et al [24,[52][53][54][55] clearly has some parallels with Gradient-Descent or Nelder-Mead search techniques. However, one can also see that this correlation can only be a part of the whole story and that GAs must confer their own advantages: indeed a very strong correlation would render the genotype representation superfluous because the dynamics of the flow could be mapped directly to the phenotype. Some power of GAs clearly lies then in the fact that small changes in the genotype can result in huge changes in the phenotype, so that the GA samples the whole space before it begins to converge. Holland's schema theory idealised this behaviour in terms of hidden building blocks within the genotype. A schema is a representation of some crucial set of digits that confers some favourable characteristic, and in this instance might look something like S = 1 * * * 0 0 * . This example has 3 entries that we are interested in (hence we say it is order 3) and 4 entries that we do not care about, which are labelled with a wildcard * . The probabilistic arguments of Holland (summarised in [7]) then suggest that schemata of this kind are important because selection favours the propagation of shorter strings of data: small subsections of the genome that confer fitness dominate first and, once they are shared by the majority of the population, crossover does not affect them. As a result, so the theory goes, if the schema S confers fitness on an individual, then it will grow exponentially within the population until an equilibrium is reached where the attraction to the solution balances against the repulsion due to the mutation.
Thus the overall behaviour suggested by schema theory is that a GA incorporates and balances competing forces. Selection and breeding tends to produce convergence around local maxima in the fitness landscape, drawing the population in over generations. On the other hand the effect of mutation is to push the population away from local maxima (on average), so that as a whole it can explore the entire parameter space. The power of GAs is then that they are sensitive to the entire landscape, but simultaneously can respond to and converge on interesting regions.
The schema argument fell temporarily out of favour in the 1990s, for reasons which are summarised in Ref. [30], mainly on the grounds that it does not consider the constructive effects of cross-over, and that it assumes a constant contribution to the fitness of a given schema. In reality, the population is not static, so that changes in the rest of the genotype tend to alter the relative fitness of the schema after the first few generations in a run. Despite this objection, the idea was rehabilitated by Poli et al in a suite of papers (see Ref. [30] for a comprehensive review) which developed so-called exact schema based models, and showed their relation to Markov chains. Simultaneously there have been several other approaches to understanding GAs, for example statistical mechanical approaches [58] or stochastic differential equations [59], which each seem to have their own drawback. In the former approach, one assigns the fitness function with a Boltzmann weighting, while the latter approach is restricted to simple cases such as the so-called onemax problem which is equivalent to solving the one-dimensional Ising model. In these cases it is likely that the set-up of the configuration influences the results so that the conclusions are not general but really only applicable to the very specific configurations used to frame the argument. Despite these limitations, the schema theory remains as one of the best interpretations of the behaviour of the GA.

Heterotic models with monad bundles
Having presented the main features of GAs, we now turn to a concise discussion of the geometric ingredients entering the construction of E 8 × E 8 heterotic string compactifications on Calabi-Yau three-folds with monad bundles. More details can, for example, be found in Refs. [10,[60][61][62][63][64][65]. In particular, we will specify the degree of similarity to the properties of the standard model, that qualifies a particular heterotic model to be designated as "perfect". From the perspective of the GA this amounts to working out the phenotype.

Heterotic CY models
As already alluded to, heterotic E 8 ×E 8 compactifications are based on the data (X, V,Ṽ , C), where X is a Calabi-Yau three-fold, V andṼ are slope poly-stable vector bundles over X whose structure groups embed into E 8 and C is an effective curve class on X. Poly-stability is the mathematical condition to ensure the bundles preserve supersymmetry and it is, in general, difficult and timeconsuming to check explicitly. For the purpose of our algorithmic approach we will only check a necessary condition known as Hoppe's criterion. It says that a poly-stable bundle V over X must necessarily satisfy Another constraint on the bundle V arises from the fact that only special unitary groups can be embedded into E 8 . This means that we require Anomaly cancellation: The sectors of the four-dimensional theory which originate from the bundleṼ and the five branes wrapping a curve with class C are hidden sectors, in the sense that they couple to the observable sector associated to V only gravitationally. Since we are interested in the observable particle physics, these sectors will not be constructed explicitly. However, for a consistent model, we need to satisfy the anomaly cancellation condition which relates the various sectors of the theory. This is done by demanding that which guarantees the existence of a supersymmetric completion of the model with a trivial hidden bundleṼ and five-branes wrapping a curve with class C = c 2 (T X) + ch 2 (V ).
Spectrum: For the purpose of this paper, we will consider "observable" bundles V with structure group SU (4) which breaks one of the E 8 symmetries to an SO(10) GUT group. The SO(10) multiplets arise from the decomposition of the 248 adjoint of E 8 under the maximal subgroup SU (4) × SO (10), which reads gauge families anti-families Higgs bundle bosons moduli The number of each multiplet is given by the first cohomology of V , V * and its various products, as indicated in the last row above. Of particular importance is the chiral asymmetry of families which can be expressed in terms of the index of V as The last equality in the above formula assumes that h 0 (V ) = h 3 (V ) = 0 which follows from the poly-stability of V , via Hoppe's criterion (3.1).
Equivariance: The SO(10) GUT symmetry needs to be broken to the standard model group and this is done by including a discrete Wilson line into the construction, which is possible when the manifold is non-simply connected. However, most of the typical constructions lead to Calabi-Yau three-folds that are simply connected. On the other hand, if a freely-acting symmetry Γ on X exists, this can be used to produce a quotientX = X/Γ with non-trivial fundamental group. The Wilson line is then included on the resulting quotient manifoldX = X/Γ. In order to end up with three chiral families after performing this quotient, the upstairs model on X requires a chiral asymmetry of where |Γ| is the order of the discrete group Γ. For this construction to make sense, the upstairs bundle V → X needs to descend to a bundleV →X on the quotient CY and this is the case iff V admits a Γ-equivariant structure. This is typically a non-trivial constraint which will be discussed in more detail below, when we introduce monad bundles.

Calabi-Yau three-folds
To realise the above construction we require, in a first instance, explicit CY manifolds. It turns out that, for group-theoretical reasons, the Wilson-line breaking of SO(10) to the Standard Model requires a discrete group Γ which is at least Z 3 × Z 3 . While the number of known CY manifolds is large, only a handful of known examples have symmetry groups this large [10,[66][67][68], so the choice of manifold is, in fact, rather constrained 1 . For the purpose of this paper, we will focus on two of these manifolds, both realised as complete intersections in products of projective spaces (CICYs) and represented by the following configuration matrices These manifolds are also known as the bi-cubic CY and the triple tri-linear CY, respectively. They are embedded in the ambient spaces A = (P 2 ) ×2 for the bi-cubic and A = (P 2 ) ×3 for the triple tri-linear. The manifold itself is cut out by the common zero locus of homogeneous polynomials whose multi-degrees correspond to the columns of the above matrices. The superscripts provide the Hodge numbers h 1,1 (X), h 2,1 (X) and the subscript the Euler number, η(X), of the manifold. Both manifolds have a freely-acting symmetry Γ = Z 3 × Z 3 , as required for the SO(10) GUT breaking. We will also need the second Chern class c 2i (T X) of the tangent bundles and the triple-intersection numbers d ijk . Relative to the obvious basis (J 1 , . . . , J h ) of H 2 (X), obtained by restricting the Kähler forms of the projective ambient space factors to X, and the associated dual basis of H 4 (X), they are given by (3.8)

Monad bundles
Explicit bundle construction frequently rely on holomorphic line bundles as basic building blocks. Holomorphic line bundles are classified by their first Chern class. Relative to an integral basis (J 1 , . . . , J h ) of H 2 (X), the first Chern class of a line bundle L → X can be represented by an integer vector k ∈ Z h , such that c 1 (L) = k i J i , hence the notation L = O X (k).
Monad bundles are defined by a short exact sequence where, for our purposes, B and C are each taken to be line bundle sums. Thus such monad bundles are specified by a pair (B, C) of line bundle sums which can be identified with an integer matrix (b 1 , . . . , b r B , c 1 , . . . , c r C ) of size h × (r B + r C ). The bundle V is isomorphic to the kernel, Ker(f ) of the monad map f which can be thought of as a r C × r B matrix with polynomial entries A given choice of (B, C) might lead to a monad sheaf, rather than a monad bundle. This happens when the degeneracy locus of f , that is, the locus on X where the rank of f is less than maximal, is non-trivial. Since we are considering smooth, geometrical models we should avoid sheafs, so we have to check "bundleness" of V by computing the dimension Typically, the matrix f exhibits a pattern of zero and non-zero entries, depending on whether Γ(O X (c a − b i )) is trivial or non-trivial, and for otherwise generic choices, the dimension d deg can be determined from this pattern. The model is acceptable iff d deg = −1.
A monad map with too many zero entries might also lead to a split bundle, that is, to a bundle V whose structure group is a proper sub-group of SU (4), for example S(U (2) × U (2)). These split bundles lead to a GUT group larger than SO(10) and we do not consider this possibility here. To this end, we introduce the number n split = number of splits of the bundle V As mentioned earlier, Γ-equivariance of V is a non-trivial constraint but a full check based on presently available methods is not compatible with an efficient computational realisation. Instead, we will perform a strong necessary check for B and C to possess a Γ-equivariant structure. This is done by checking that each unique line bundle L in B or C has an index divisible by |Γ|. If a line bundle L appears with multiplicity m in either B or C, we require that ind(L ⊕m ) is divisible by |Γ| . To check the other constrains we require the following formulae for the Chern classes of a monad bundle: (3.14) The conditions appearing on the right-hand sides of these equations originate from the requirement that V has an SU (4) structure group, from the anomaly cancellation condition (3.4) and from the constraint (3.6) on the chiral asymmetry.
Calculating cohomology for monad bundles requires the long exact sequence in cohomology associated to the monad sequence (3.9). Such computations can require ranks of maps and can be time-consuming. For this reason, checking Hoppe's criterion for poly-stability in full is not feasible while running a GA. Instead, we will use that Hoppe's criterion is violated and, hence, that V is non-supersymmetric if Fortunately, analytic formulae for zeroth line bundle cohomology exist for both manifolds under considerations [10], so the above expressions can be checked efficiently.

Monad bundles and genetic algorithms 4.1 The environment
Let us now outline the "environment" that will form the basis of our search. We will consider models consisting of monad bundles (3.9), specified by matrices (B, C), on a fixed CY three-fold X, in practice one of the two manifolds in Eq. (3.7). The integer entries in (B, C) are not, a priori, bounded but will be restricted as so that the environment becomes finite. Given that we are aiming to match relatively small integers, such as the required index, this does not seem to be a serious limitation. In fact, model building experience suggests [6,10,69,70] that viable models typically arise for relatively small integer entries. The first Chern class constraint, c 1 (V ) = 0, can be easily solved (see Eq. (3.14)), for example by fixing the last line bundle in C in terms of the others. For this reason, we constrain the environment to pairs (B, C) that satisfy c 1 (B) = c 1 (C). If we allowed, say, 10 values per entry, the number of states in the environment is of order 10 h(r B +r C −1) .

(4.2)
Since r B + r C − 1 ≥ 5 and we consider h = h 1,1 (X) ∈ {2, 3} this is quite sizeable and not suited for systematic scanning. Model building experience also shows [6,10,69,70] that the fraction of viable models within this environment is so small that a purely random search could not be successful.
The environment has a large degeneracy, governed by the symmetry where the first factor is H = S 2 for the bi-cubic CY and H = S 3 for the triple tri-linear CY and arises due to the symmetries of the configuration matrices; it corresponds to permutations of the rows of the matrix (B, C). The permutation groups S r B and S r C arise from permutations of the line bundles (columns) in B and C which of course do not affect the resulting monad bundle V . For the purpose of applying GAs, we have kept this degeneracy, as the ordering required for its elimination would significantly complicate the algorithm. Indeed the very notion of ordering by genotype, which would have to be performed even after each mutation, is somewhat at odds with the principle of GAs although it could in principle be included. However, redundancy in the models will be ultimately be eliminated from the list of perfect models found by the algorithm.    Table 1, such a state allows for a supersymmetric, anomaly-free completion, it is based on an SO(10) GUT symmetry, has the correct chiral asymmetry of families and it passes the checks for bundle equivariance and poly-stability discussed earlier. Of course such perfect models might then still fail at a fairly elementary level, for example due to the presence of undesirable 16-multiplets, a failure of stability or a failure of equivariance. Checking these properties requires cohomology calculations which, with the available methods in commutative algebra, are too time-consuming to be carried out during the GA evolution. However, experience from "by-hand" model building [6,69,70] and from RL [10] shows that a significant fraction of the perfect models also satisfy the more stringent cohomological constraints. In other words, while some of the conditions in Table 1 are only necessary, they are still sufficiently strong 2 .
The monad environment we have just described was also used in the context of RL [10], where the value function in Table 1 was used to determine the reward. It was shown in Ref. [10] that RL is extremely successful on this environment, with trained policy networks efficiently leading to perfect models for 100% of episodes, and many, previously unknown, models being discovered. Using the same environment allows a direct comparison between RL and GAs, which we carry out in the next section. In summary then, our environment is the same as that used for RL in Ref. [10], which has been realised as a MATHEMATICA [71] package.

Genetic algorithm
To realise the GA we choose the range of entries such that b max − b min = 2 n B − 1 and c max − c min = 2 n C − 1, representing every integer in B and C in a binary encoding, with n B and, respectively, n C bits, resulting in a bit sequence of total length = h(r B n B + (r C − 1)n C ). The typical size of a population consists of N pop 100 − 300 states, each represented by such a binary encoding of the matrix (B, C). The population is then evolved by crossing and mutation, as discussed in Section 2.
With the fitness function specified by the intrinsic value v(B, C), we found that rank-weighting selection, elitism and a simple one-point cross breeding work well. The GA was implemented in two ways, as a MATHEMATICA [71] and as a Python package, in both cases being coupled to the MATHEMATICA package realising the environment. The two realisations lead to similar results. Working with integers that can be encoded on 3 bits is convenient for the purpose of comparing the GA and RL techniques. However, it is possible and in fact easy to further increase the range of line bundle integers, as the corresponding augmentation in the computational time is not significant.
As mentioned above, the environment has a large degeneracy. Equivalent models arise from permuting the two P 2 -factors in the bicubic embedding, as well as from permuting the line bundles in B and C. This amounts to a group of order 2! · 6! · 2! = 2800. Of course, this does not necessarily imply that every perfect model can be found in 2800 different sites, since some permutations may leave particular matrices (B, C) unchanged, but it nevertheless gives some idea of the degree of degeneracy present in the environment.
It is useful to illustrate the performance of the GA graphically. We set the size of the population to N pop = 250 randomly initialised states and let the GA run for 200 generations with a mutation rate of 0.004 and the parameter α from Eq. (2.3) set to α = 3. The time taken by such a run is relatively small (<100 seconds on a laptop). Figure 1 shows a typical evolution of the population.
The histogram on the left shows how the population progresses towards greater fitness (intrinsic value). The plot on the right shows the fraction of perfect models in the population, that is, the fraction of models with a value v(B, C) = 0. By generation 150 more than half of the population corresponds to perfect models. For this run, a total of 200 × 250 states have been visited, many of them multiple times. Of these 12, 665 correspond to perfect models with 48 being distinct. After eliminating redundancies, 18 non-equivalent models remain.
The fraction of perfect models can be easily increased by lowering the mutation rate. However, as one might anticipate from the discussion in the introduction, this also induces a degree of stagnation, with a tendency for the same states to emerge repeatedly.
A few remarks are in order. Firstly, to appreciate the performance of the algorithm, note that the  By comparison, the number of states visited in the above run, namely 50, 000, represents only a tiny fraction of the space. However, the GA was capable of finding 48 perfect states, while a random search over millions of states would typically lead to no perfect states at all.
Secondly, as already noted, the GA has a tendency of visiting the same states multiple times. It is interesting to plot the total number of perfect states found after n generations as a function of n. For our illustrative run such a plot is shown in Figure 2, which suggests that there is no additional benefit in letting the population evolve beyond a certain generation (n 150). Finally, it is useful to compute the degeneracy of the 18 states that remain after removing redundancies. By performing all the allowed permutations a number of 19, 080 states are obtained. This should be compared with the product 18 × 2800 = 50, 400 which turns out to be an overestimate by more than a factor of 2. Moreover, what this computation indicates is that a single run of the GA is not enough if the aim was to find all the perfect states available in the environment and that 1, 000 further runs, which would take about 1 day on a single machine, would be just about enough to find the other redundant representations of the 18 states found in the first run. Of course, many more new states, not related to the 18 by permutations, would be expected to arise in such a search, which implies that a comprehensive search would require several, possibly tens of thousands of GA runs. With 10, 000 runs this would amount to exploring ∼ 0.01% of the environment.

A quasi-comprehensive GA search
The illustrative results presented above show that GAs can be efficiently used to identify monad bundle models on the bicubic manifold that satisfy the phenomenological requirements discussed in Section 3. A natural question to ask is whether the GA can find all (or almost all) perfect models available within the environment. Answering this question is difficult, mainly because the environment is too large to be subjected to a systematic scan which would ensure that no perfect state was missed. The computational time needed for a systematic scan within the environment bounds can be easily estimated to about 200 core years, which although attainable, is by no means something that one would easily invest in.
A partial answer to the question can be achieved by making use of the redundancy present in the environment. The procedure is as follows. We perform 10,000 runs keeping the same settings as above, namely N pop = 250, the length of a run equal to 200 generations, a mutation rate of 0.004 and α = 3. The required computation time for this is about 10 core days 3 . We extract the number of perfect states found after n runs and then ask how many of these are left after removing redundancies. Now, if the GA is capable of finding almost all the states, then the expectation would be that the redundancy reduced number of perfect states after n runs would saturate as a function of the total number of perfect states found after n runs. In other words the GA is only finding versions of models that it has found before. The plot in Figure 3 confirms this expectation.  The search lead to a number of 96, 705 perfect states, out of which 639 are inequivalent (not related by a permutation). Performing all possible permutations on the 639 states leads to 481, 680 perfect states that must be present in the environment. If allowed more time, the GA would of course be expected to find all of these. However, the saturation plot in Figure 3 as well as the plots in Figure 4 suggest that there are not many more inequivalent states to be found. We take this as an indication that the search has already successfully found most of the inequivalent perfect states.

Enlarging the search space
So far the minimal and the maximal values of the entries in B and C were fixed such that all integers can be encoded on 3 bits. Upgrading to integers that can be encoded on 4 bits leaves the  To answer this question we need to see if the perfect models found are well within the boundary of the environment or if a significant fraction of them lies close to the boundary. The bar charts in Figure 5 show that in fact most of the models lie on the boundary (for instance, 451 out of the 639 inequivalent models include the integer 4 in B). This suggests the need for increasing further the entries in B and C. However, before embarking on such an enterprise, it is worth looking deeper into the quality of the monad bundles found. The requirement of stability imposes that h 0 (V ) = h 3 (V ) = 0. Moreover, in order to avoid the presence of 16-multiplets, the bundle has to satisfy h 2 (V ) = 0. These cohomological constraints are too computationally intensive to have included them in the search, however they can be carried out on the set of 639 inequivalent models found on the bicubic manifold.
Moreover, many of the models found contain a repeated line bundle in B and C. While this is not a problem in itself, such models are equivalent to monad bundles with rk(B) = 5 and rk(C) = 1. We eliminate these from the discussion, which leaves 214 models of genuine rk(B) = 6 and rk(C) = 2 type. Out of these, 67 have the right cohomologial properties mentioned above. Further stability models found (inequivalent) models of genuine type (6,2) models of genuine type (6,2) satisfying cohomological constraints stable models of genuine type (6,2) satisfying cohomological constraints checks reduce this number to 29 models. The distribution of entries in B for the remaining 29 models is shown in blue in Figure 6 and indicates that most of the models have all the line bundle entries in the range {−2, −1, 0, 1, 2}. In view of this, the expectation is that very few, if any, additional models satisfying the additional constraints would be found by further enlarging the size of the environment.

GA scan versus RL scan
In Ref. [10] a similar search for physically viable models on the bicubic manifold with monad bundles was performed using methods of reinforcement learning (RL). The same fitness function v(B, C) as described above was used in that context to compute the intrinsic value of states (B, C). When moving from one state to another, the RL agent receives a reward equal to the difference between the intrinsic values of the two states if the difference was positive or a fixed penalty if the difference was negative. The search is divided into episodes of fixed maximal length and along each such episode the return G t of a state s t is computed as where r t is the reward (penalty) associated with the step s t → s t+1 and the discount factor γ is a subunitary number close to 1. The aim of the game is then to devise a strategy of finding a path that would maximise the return of any randomly chosen starting point, by eventually landing on a perfect state (also called a terminal state in the RL context).
More details can be found in Ref. [10]. What is important for the present discussion is that, although the RL and GA methods rely on the same intrinsic value function, their underlying philosophies are very different, which makes for an interesting comparison.
For the RL we use a policy network trained on episodes of maximal length 64. After 100, 000 rounds of training, a phase which takes several hours on a laptop, the RL agent is capable of finding terminal states from virtually any random starting point within episodes of average length ∼30.
The training curves are shown in Figure 7. In order to get a preliminary assessment of the performance of the trained network we let it run for 10, 000 episodes which takes about 16 minutes on a single machine. The run produced 2460 terminal states, out of which only 72 are distinct. The repetitions are not evenly distributed over the 72 distinct states. In fact, one of these states accounts for almost half of the terminal states found. Incidentally, this particular state corresponds to a monad bundle that can be reduced to a rk(B) = 5 and rk(C) = 1 case. The genuine rk(B) = 6 and rk(C) = 2 are found with much smaller frequencies, most of them appearing only once in the dataset. Thus, although the network is capable of finding terminal states, a small number of these sates have very large basins of attraction and are found much more often than the others. This is a source of inefficiency for the RL implementation.  We let the trained network search for terminal states from 1.7·10 9 random initial states. This search, which took 35 core days produced ∼10, 000 perfect states, out of which 643 were inequivalent. By comparison, the computational time for the GA scan was 10 core days and lead to ∼100, 000 perfect states, out of which 639 were inequivalent. The saturation curve of the number of inequivalent models as a function of the total number of models found at various stages of the RL search is shown in Figure 8 together with the corresponding curve for the first 10, 000 perfect states found during the GA search. Figure 9 shows the progression of the number of (inequivalent) models found with the number of states visited.  A number of remarks can now be made. Firstly, the GA search is more efficient at finding perfect states than the RL search. Roughly the same number of inequivalent perfect states have been found in the two searches, though the RL search took 3.5 times longer than the GA. Secondly, the GA is more prolific, by more than an order of magnitude compared to the RL search, at finding redundant perfect states. This is very likely a by-product of crossover, which is prone to producing redundant states by permuting various parts of the genome. Note that various techniques such as crowding penalties might be considered to remove this inefficiency.
Finally, there is a substantial overlap between the inequivalent models found in the two searches, with only ∼50 models in each complement. This brings the total number of inequivalent models to 689. The fact that two such different methods lead to virtually the same dataset of models is another even stronger indication that it has been possible to achieve a high degree of comprehensiveness by scanning only a tiny fraction of the environment.
While both methods produce nearly identical data sets of perfect models after sufficient running time their trajectories during evolution/training appear to be quite different. Using nonlinear visualisation techniques such as Sammon maps we find that the set of perfect models found at an early stage is quite different for the two methods. In particular, as illustrated by Fig. 8, GAs tend to produce a significantly larger degree of redundancy than RL.
6 Monad bundles on the triple trilinear manifold 6.1 A typical GA run on the triple trilinear As in the previous section we consider monad bundles with rk(B) = 6 and rk(C) = 2 and set b min = −3, c min = 0 and b max = 4, c max = 7, such that each integer entry in the matrix (B, C) can be represented by a 3-bit binary sequence. For the triple trilinear manifold the degeneracy of the environment is even larger than in the case of the bicubic, due to the six possible permutations of the three P 2 -factors in the ambient space. The corresponding group has order 8640.
As in the previous case, we illustrate the performance of the GA graphically. We keep the size of the population to N pop = 250 randomly initialised states. The GA needs to run a little longer until a significant fraction of the population corresponds to perfect states. Evolutionary episodes of 400 generations with a mutation rate of 0.004 and α = 4 seem to work well. The time taken by such a run is relatively small (∼450 seconds on a single machine).   Figure 10 shows a typical evolution of a population of models on the triple trilinear manifold. By generation 300 about half of the population corresponds to perfect models. For this run, a total of 400 × 250 states have been visited, many of them multiple times. Of these 15, 377 correspond to perfect models, with 43 being distinct. After eliminating redundancies, 41 non-equivalent models remain.
The size of the search space is 8 21 10 19 .
which is orders of magnitude larger than the bi-cubic space. The number of states visited during the above run, namely 100, 000 states (not necessarily distinct), represents a minuscule fraction of the space. Nevertheless, it is remarkable that 43 perfect states were found, and these perfect states lead to 41 inequivalent models. The total number of perfect states found after n generations as a function of n is plotted in Figure 11, which suggests that the number of generations, 400, is suitably chosen given the other hyper-parameters. This continued efficacy despite the much larger search space suggests that the difficulty of finding a solution is increasing only polynomially for the GA as one might have hoped.

GA search on the triple trilinear
The environment being much larger than in the case of the bicubic manifold, we have not attempted to reach the same degree of comprehensiveness as before. After running for 15 core days, the GA found ∼ 100, 000 perfect states, three quarters of which are inequivalent. Figure 12 shows that the number of inequivalent states is far from saturation at the end of the run, suggesting that many more other inequivalent perfect states are still to be found within the environment.

Conclusion
In this paper we have explored the efficacy of genetic algorithms (GAs) as a search method for probing the landscape of four-dimensional N = 1 heterotic compactifications. GAs turn out to be very successful for engineering the compactification data and for finding perfect models, that is, models which resemble the Standard Model. The spaces of models scrutinised here were extremely large (∼ 10 12 in the case of the bicubic manifold and ∼ 10 19 for the triple trilinear manifold). In the case of the bicubic manifold a systematic scan, although feasible, would be extremely time consuming, while for the triple trilinear manifold, a systematic scan is practically unfeasible. Despite this, the GA was able to find tens of models within the first few minutes of running on a single machine. The method is not only efficient, but also highly parallelisable.
After 10 days of running, the GA accomplished a high degree of comprehensiveness in finding all the perfect models available on the bicubic manifold. Indirect evidence in support of this claim was obtained by exploiting the degeneracy of the environment, as well as by comparison of the results of the GA search with those obtained through reinforcement learning (RL). Although the two methods are very different, they lead to datasets of inequivalent models that have an overlap of over 90%. This suggests that the details of the optimisation process are not really essential once the processes begin to saturate, provided that they share the same incentives.
Our comparative analysis of the two search methods revealed that GA is, overall, more efficient by about an order of magnitude in finding perfect states than RL. However, GA is also more prone to finding models that are related by redundant permutations, a feature which is likely related to the crossover mechanism. In principle techniques such as crowding penalties (with the redundancy incorporated) might improve this aspect. An important factor contributing to the lower output of the RL search was the fact that certain perfect states had very large basins of attraction and consequently were found many times, while other perfect states had small basins of attraction and needed many trial episodes to be found.
The combined RL and GA searches on the bicubic manifold led to some 700 inequivalent perfect models. Only a tenth of these correspond to monad bundles of genuine (rk(B), rk(C)) = (6, 2) type with the correct cohomology dimensions. Further stability checks reduce the number of viable models by another factor of 2 -the corresponding monad bundles have been tabulated in Appendix A. Unfortunately these checks had to be performed at the end of the search due to the limitations imposed by cohomology computations. However, it may be the case that a quick method for computing monad bundle cohomology can be developed, for instance in the form of an analytic formula similar to those found for line bundle sums. If possible, this would facilitate the implementation of further checks at the search stage, which would significantly improve the performance of both GAs and RL.
Improving the performance of the search methods is of particular importance for large environments, such as that corresponding to the triple trilinear manifold. In the absence of more stringent notion of perfect model, achieving any degree of comprehensiveness is difficult in such large environments; consequently, partial searches such as the one presented in Section 6.2 may miss many relevant models.
Looking more broadly, one is irresistably drawn to the question of whether the optimisation processes that we have been exploring here might have a bearing on physics in the early Universe, in the sense that a qualitatively similar kind of selection may have happened involving, for instance, topological and small-instanton transitions between different vacua. Of course it is not known what (if any) incentives might favour the physics that we observe in our world, although there have over the years been varyingly speculative suggestions (and we do not need to add to the list). But if this is the case and if the idea of a multiverse is correct, then the results that we observe with GA and RL searches do at least support the idea that selection mechanisms can successfully operate within the vast search space, and indeed that they do so in a quite convergent way. This suggests the possibility of alternative pictures, quite distinct from that of a randomly populated multiverse, in which the properties of this world would be more common than one would naively expect.

A Monad bundles on the bicubic
The table below contains the 28 monad bundles of genuine type (r B , r C ) = (6, 2) obtained on the bicubic manifold X, leading to models with 27 families and no anti-families. For these bundles equivariance with respect to the Z 3 × Z 3 action on the manifold has been checked at the level of index divisibility for the line bundles in B and C. A number of necessary stability checks have been carried out. More precisely, we have found all the line bundles with entries in the range [−3, 3] that inject into V and have checked that there exists a non-empty locus in the Kähler cone where all these line bundles have a negative slope, thus leaving V slope poly-stable.
In order to save space we use the notation O(k 1 , k 2 ) instead of O X (k 1 , k 2 ). All the line bundles are to be understood on the bicubic threefold.