Combining Evolutionary Strategies and Novelty Detection to go Beyond the Alignment Limit of the $Z_3$ 3HDM

We present a novel Artificial Intelligence approach for Beyond the Standard Model parameter space scans by augmenting an Evolutionary Strategy with Novelty Detection. Our approach leverages the power of Evolutionary Strategies, previously shown to quickly converge to the valid regions of the parameter space, with a \emph{novelty reward} to continue exploration once converged. Taking the $Z_3$ 3HDM as our Physics case, we show how our methodology allows us to quickly explore highly constrained multidimensional parameter spaces, providing up to eight orders of magnitude higher sampling efficiency when compared with pure random sampling and up to four orders of magnitude when compared to random sampling around the alignment limit. In turn, this enables us to explore regions of the parameter space that have been hitherto overlooked, leading to the possibility of novel phenomenological realisations of the $Z_3$ 3HDM that had not been considered before.


I. INTRODUCTION
The Standard Model (SM) of particle physics has demonstrated remarkable success in accurately describing the electroweak and strong interactions.However, experimental challenges such as neutrino mass, dark matter, and the baryonic asymmetry of the Universe prompt the exploration of physics Beyond the SM (BSM).In many cases, BSM theories involve expanding the minimal scalar sector of the SM, characterised by a single Higgs doublet.
Multi-Higgs doublet models are particularly prominent among these extensions, primarily because they maintain the tree-level value of the electroweak ρ-parameter, in good agreement with experimental observations.The extensively studied two Higgs doublet model (2HDM) [1] has provided valuable insights.More recently, there has been a surge of interest in the investigation of three Higgs doublet models (3HDMs) [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], where the scalar sector encompasses three Higgs doublets.These models show great promise as they possess the essential ingredients to address the challenges posed by dark matter and the baryonic asymmetry of the Universe.
With no unambiguous signs of New Physics in general and of extra exotic scalars in particular, BSM phenomenology is faced with an ever increasing and ever more restricting list of (direct and indirect) experimental constraints, in addition to any theoretical and selfconsistency constraints.From a model building perspective, this means that the allowed regions of BSM parameter spaces are effectively shrinking, making finding such regions ever more difficult or outright impractical, even if those regions host points which are very good fits to the data.To mitigate this, BSM phenomenologists often simplify the problem in two possible ways.The first approach is to simplify the constraints and leave some questions unanswered, to be eventually addressed by an Ultra Violet completion.The second approach relies on simplifying the sampling space by changing the priors from which the parameters are drawn, usually by restricting the parameter space to a subregion where known valid points had been found before or for which there is a limiting case, such as the case of the alignment limit in multi Higgses models.While in the former case one ends up with an incomplete model, in the latter case one is left with the worrisome prospect of missing out possible phenomenological signatures.
In recent years, Artificial Intelligence (AI) in general, and Machine Learning (ML) in particular, have received considerable attention by the High Energy Physics (HEP) community with a wide range of applications [17].Of particular interest to this work are the ongoing attempts to explore the parameter space of highly constrained and multidimensional BSM scenarios, where sampling points that respect theoretical and experimental bounds poses a great challenge.The first attempts to mitigate this problem using AI/ML are based on using supervised classifiers and regressors to produce estimates of BSM predictions of physical observables and physical quantities [18][19][20][21][22], bypassing the computational overhead often associated with these.Other approaches leverage the active learning framework to track the ground truth of the observables to draw the boundary of the allowed regions of the parameter space [23,24].However, we point out that the usage of AI/ML for BSM parameter spaces studies is not restricted to exploration, as generative models have been studied to provide a possible way of replicating points from valid regions [25] or to gain new insights through latent space visualisation [26].
The approaches cited above rely heavily on the quality of the trained ML model, namely on the coverage of the parameter space provided by the points used, especially in or near valid regions.A different approach was presented in [27], where different AI-based exploration algorithms were used to explore the parameter space of the cMSSM (four parameters) and the pMSSM (16 parameters) by reframing the problem as a black-box optimisation problem.
In such an approach, exploration starts with just a few random points from which iterative algorithms will progressively suggest better points through trial and error.Although the Physics cases in [27] were not especially realistic, as only Higgs mass and Dark Matter relic density constraints were used, the methodology provided orders of magnitude of sampling efficiency improvement over random sampling, while still capturing a global picture of the parameter space.
In this paper, we will extend and build on top of that of [27].In that work, an Evolutionary Strategy algorithm was observed to eagerly converge to the valid region of the parameter space and to stop exploring once converged.To mitigate this, in [27] the algorithm was endowed with a restart strategy to draw a more global picture of the valid region of the parameter space.In this work, we will take a different approach by incorporating a novelty reward into the black-box optimisation problem to drive the exploration algorithm away from the regions already explored.As we will see, our approach retains the benefits of drastically improving sampling efficiency while gaining a novel approach to charter the valid region of the parameter space.
Although the proposed approach is general to any BSM parameter space,1 we will use our novel methodology to perform a realistic search on the highly constrained Z 3 3HDM model, where all known and possible constraints will be considered.This poses a terrific challenge for sampling valid points from the parameter space, which has led previous studies to consider sampling around the alignment limit.Therefore, the Z 3 3HDM provides an ideal scenario not only to develop our methodology, but to use to do new Physics by exploring points beyond random sampling strategies around the alignment limit and to discover novel phenomenological realisations of the Z 3 3HDM that have been obfuscated in the past such strategies.
This paper is organised as follows.In section II we present the Z 3 3HDM model which is the Physics subject of our study.In section III we present its constraints, both theoretical and experimental.In section IV we outline the random sampling strategy near the alignment limit, which we will use to compare our methodology.In section V we introduce the AI scan strategy, redefining the scan as black-box optimisation, and the novelty award based on density estimation.In section VI we present and analyse the results obtained with our methodology, showcasing the versatility of our approach.Finally, in section VII we conclude our discussion and point out novel directions of work.

A. Scalar Sector
For the potential of the Z 3 3HDM model, denoted V Z 3 , we use the conventions of [6].The Z 3 -invariant terms, given by ϕ i → ϕ ′ i = (S Z 3 ) ij ϕ j , can be expressed as: with the quartic part represented by: which includes terms, m 2 12 , m 2 13 , and m 2 23 , responsible for breaking the symmetry softly.Following spontaneous symmetry breaking (SSB), the three doublets can be parameterised in terms of their component fields as: where v i / √ 2 corresponds to the vacuum expectation value (vev) for the neutral component of ϕ i .It is assumed that the scalar sector of the model explicitly and spontaneously conserves CP.Under this assumption, all parameters in the scalar potential are real, and the vevs v 1 , v 2 , v 3 are also real.
The scalar potential in eq.(1) contains eighteen parameters, and the vevs are parameterised as follows: leading to the Higgs basis [28][29][30] obtained by the following rotation, where here we have used the short notation, c x ≡ cos x, s x ≡ sin x.
Orthogonal matrices, denoted as R, P and Q, diagonalise the squared mass matrices in the CP-even scalar, CP-odd scalar, and charged scalar sectors.These matrices are crucial for transforming the weak basis into the physical mass basis for states with well-defined masses.Although this has already been discussed before [6,10,31], for completeness and to fix our notation, we give here the rotations that relate the mass eigenstates to the weak basis states.
For the neutral scalar sector, the mass terms can be extracted through the following rotation, where we take h 1 ≡ h 125 to the be the 125 GeV Higgs boson found at LHC.The form chosen where the matrices R i are For the charged and pseudoscalar sectors, the physical scalars can be obtained via the where, the rotation matrices are given by For later use, we define the matrices P and Q as the combinations that connect the weak basis to the physical mass basis for the CP odd and charged Higgs scalars, respectively, As the states in the physical basis have well-defined masses, we can obtain relations between the set and the potential parameters in eq.(1), as shown in Ref. [6,10,31].

B. Higgs-Fermion Yukawa Interactions
In the Type-I models considered here2 , fermion fields are unaffected by the Z 3 transformation, allowing them to couple only to ϕ 3 .The Yukawa couplings to fermions are expressed compactly as: where h j ≡ (h 1 , h 2 , h 3 , A 1 , A 2 ) j represents the physical Higgs fields.For completeness, we list the couplings a f j and b f j here [12].We have, for all leptons and down quarks, for all leptons and down quarks, For the charged Higgs, H ± 1 and H ± 2 , the Yukawa couplings to fermions are expressed as: where (ψ

III. CONSTRAINTS
In this section, we outline the various constraints necessary to impose theoretical and phenomenological consistency on the model parameters.The specifics of these constraints in the context of 3HDM are well established, as documented in previous works [10,12,16].
For brevity, we provide a brief list here, deferring further elaboration to appendix A.
From a phenomenological standpoint, our primary objective is to ensure the existence of an SM-like Higgs, identifiable as the scalar boson detected at the LHC.As demonstrated in [6], achieving this involves staying close to the 'alignment limit,' characterised in the 3HDM by the conditions: In this limit, the lightest CP-even scalar, denoted as h, exhibits exact SM-like couplings at the tree level, automatically satisfying constraints from Higgs signal strengths.However, our interest lies in exploring permissible deviations from this precise alignment.To this end, we use the signal strength formalism, comparing the results with the experimental limits [32].
Subsequently, we must address constraints stemming from electroweak precision parameters, specifically S, T , and U .We use the analytic expressions of Ref. [33], contrasting them with the fit values provided in Ref. [34].Notably, similar to the 2HDM scenario, we can bypass T -parameter constraints by imposing [14]: as we will explain in section IV.
Additionally, we incorporate constraints arising from flavour data, as detailed in appendix A. Direct searches at the LHC for the heavy non-standard scalars are also considered, employing HiggsTools-1.1.3following Ref.[35], which provides a comprehensive list of relevant experimental searches.
For theoretical constraints, we insist on the perturbativity of Yukawa couplings, perturbative unitarity, and BFB (bounded from below) conditions.The implementation details for these constraints can be found in appendix A.

IV. RANDOM SCAN STRATEGY
We developed a dedicated code specifically for the Z 3 constrained 3HDM3 , building upon our earlier codes [10,38,39].A thorough exploration of the parameter space was conducted using section II A. Our fixed inputs remained v = 246, GeV and m h1 = 125, GeV.Random values were assigned within the following ranges: where the last expression applies only to the soft masses that are not obtained as derived quantities (see Ref. [10] for the complete expressions).
However, this extensive scan exhibited low efficiency (as detailed in table III below).
Recognising the significance of alignment in 3HDM [6,9,10], where alignment is defined by the lightest Higgs scalar having Standard Model (SM) couplings, we imposed constraints to enhance efficiency.Initially, aligning α 1 with β 1 and α 2 with β 2 (eq.( 18)) did not yield enough good points.Ref. [9] introduced an additional condition : This, alongside eq. ( 18), led to a symmetric form of the quartic part of the potential [9,40].
In fact, these conditions simplified the potential to where is the the SM quartic Higgs coupling.All λ i vanish or are expressed in terms of λ SM .The validity of eq. ( 18) and eq.( 23) also implies that the soft masses can be explicitly solved, that is, they are no more independent parameters (see Ref. [10] for complete expressions).Now it should be clear why all such points are good points.Due to alignment, the LHC results on the h 125 are easily obeyed, whereas the perturbativity unitarity, STU and the other constraints are automatically obeyed.
To facilitate efficiency and consider deviations from strict alignment, two conditions were introduced [10].The first denoted "Al-1", allowed a percentage deviation: The condition "Al-2" was more stringent, combining Al-1 with six additional conditions: These conditions, especially Al-2, improved the generation of meaningful points beyond the SM, even with a departure from strict alignment.

V. ARTIFICIAL INTELLIGENCE BLACK-BOX OPTIMISER SCAN STRATEGY
To quickly explore the parameter space, we will employ the AI black box optimisation approach to parameter space sampling first presented in [27].In this approach, a constraint function, C, is introduced, where O is the value of an observable (or a constrained quantity), O LB is its lower bound (i.e. its lowest allowed value), and O U B its upper bound (i.e. its highest allowed value).
Here, O is obtained by some computational routine that maps the parameter space to physical quantities, where the details of such routine are irrelevant and can therefore be taken as a black-box.If the value of O is within its lower and upper bounds, C(O) returns 0, otherwise it returns a positive number that measures how far the value of O is from its allowed interval.O is dependent on the parameters of the model, θ = (α 1 , β 1 , . . . ) ∈ P (where P is the parameter space defined by eq. ( 22)), that is, O = O(θ), which implies that Therefore, the set of valid points, V, that satisfy a constraint can be defined in terms of θ as Since C(θ) is both vanishing and minimal in V, the same set can be defined through the optimisation statement Therefore, the task of finding the points in the parameter space that satisfy constraints is the same as finding the points that minimise C(O).When faced with multiple constraints on multiple observables or constrained quantities, O i , we can then combine them into a single loss function, L, which we wish to minimise where the sum runs over all the N C constraints discussed in section III, and L ≥ 0 ∀ θ with L = 0 if and only if all constraints are met.We notice that the quantity O i does not need to be observable, such as the µ if signal strengths measured by the LHC.For example, the theoretical constraints related to BFB and unitary conditions are included as O i with the respective required bounds.The ability to mix measurements, limits, and theoretical constraints in the same loss function is a key strength of this methodology.Although O i included observables and other constrained quantities, we will often abuse terminology and refer to all O i as observables and the space they span as observable space, O.
A. The Search Algorithm: Covariant Matrix Approximation Evolutionary Strategy Having framed parameter space sampling as a black-box optimisation problem, we need to choose which AI black-box optimisation algorithm to perform the task.In [27], three different algorithms were considered: a Bayesian optimisation algorithm, a genetic algorithm, and an evolutionary strategy algorithm.Each realised different balances of the so-called exploration (how much of the parameter space explored) vs. exploitation (how fast it converged to V) trade-off.In this work, we will use the evolutionary strategy algorithm, which provides the fastest convergence speed.
The evolutionary strategy algorithm in question is the Covariant Matrix Approximation Evolutionary Strategy (CMAES) [41,42].Evolutionary Strategies (ES) are powerful numerical optimisation algorithms from the broader field of Evolutionary Computing (EC).
EC algorithms are characterised by an iterative process in which candidate solutions for a problem are tested and the best ones are used to generate new solutions.In our case, the candidate solutions are points in the parameter space, and their suitability (i.e.their fitness) is measured by the loss function function eq.(31).As opposed to Genetic Algorithms, ES do not make use of genes to generate new candidate points.Instead, in ES, new candidates are sampled from a distribution, the parameters of which are set by the best candidates from previous iterations, called generations.
In CMAES, the distribution is a highly localised multivariate normal.This normal distribution is initialised with the centre at a random point in the parameter space, and its covariant matrix is set to the identity multiplied to an overall scaling constant σ.A generation of λ candidates is sampled from the multivariate normal and their fitness is evaluated by eq. ( 31).Next, the λ candidates are sorted from best to worst, and the µ best candidates are used to compute a new mean of the normal distribution, as well as to approximate its covariant matrix.Intuitively, the change in mean progresses the algorithm in the direction of steepest descent of the loss function, just like a first-order optimisation method would, and the covariant matrix approximates the (local) Hessian of the loss function, just like a second-order optimisation method would.The difference, however, is that CMAES does not compute derivatives of the loss function, and therefore it is suitable for nonconvex, illconditioned, and even non-continuous loss functions.This feature makes CMAES converge very quickly on a wide variety of optimisation problems.We warn, however, that the intuitive description of CMAES presented above hides many of its inner working details, which are not relevant for the study at hand, and we point the interested reader to the original references provided.

B. The Novelty Reward: Histogram Based Outlier System
In [27] CMAES was shown to converge very fast compared to other AI black-box search algorithms.However, it was also observed that CMAES has limited exploration capacity due to the highly localised nature of the multivariate normal from which new candidate solutions are drawn.To mitigate this, in [27] many independent runs of CMAES with restarting heuristics were performed to draw a more global picture of the valid region of the parameter space.In this work, we present a novel approach to promote exploration by adding a novelty reward to the loss function.To achieve this, we will compute the density of valid points already found and add it to the loss function as penalty.In this way, the loss function will be minimal not only when the constraints are satisfied, but also away valid regions that have already been found, which pushes CMAES to explore new regions and effectively working as a novelty reward.
The addition of a penalty to the loss function in eq. ( 31) might produce new local minima.
For example, consider the addition of various penalties, p j , each taking values This will create a competition between penalties p j , and observable optimisation C(O i ), when i C(O i ) ≃ j p j , producing local minima away from i C(O i ) = 0 and therefore spoiling our attempt to find good points.In order to prevent this, we artificially raise the value of the loss function outside the valid region by one so that such competition never arises, i.e. we consider a new version of L, L, which guarantees that the total loss function, L T , including N p penalties is still positive semi-definite, and such that for a valid point we have i C(O i ) = 0 ⇒ 0 ≤ L T ≤ 1 in a way that the density penalty does not compete with the constraints since, i.e.
for invalid points we always have L T > 1.
Having defined how penalties can be added to the loss function without spoiling our implementation of a black-box optimisation algorithm to find valid points, we now have to choose how to compute and quantify the penalty to produce the novelty reward.The first thing we need to address is how to compute the point density.This task is known in the Machine Learning (ML) literature as density estimation, and for large multidimensional datasets it is a very challenging task.Furthermore, not only we want to be able to estimate the point density accurately, we do not want the density estimation to be prohibitively slow to compute.After some preliminary exploratory experimentation, we identified the Histogram Based Outlier System (HBOS) [43] as a fast and easy to implement solution. 4BOS has also been previously explored in the context of model-independent searches for new physics using anomaly detection [44].HBOS works by fitting a histogram with a fixed number of bins, N bins , to each dimension, i.e for each parameter.A density penalty for a point θ is obtained by summing the heights, h j , of the bins on which the values of the parameters, θ j , belong. 5The penalties are normalised to be p ∈ [0, 1], so that a novelty point has penalty 0 and a point too similar to the ones already seen has maximal penalty 1.Furthermore, we notice that the penalty over the parameter space density needs not to be over all parameters, θ = (α 1 , β 2 , . . .), but can be focused on only a subset of these, {θ j } -this is especially useful to promote focused exploration in parameters of interest.
Whilst the discussion above focusses on parameter space density penalty, it can be extended to other spaces.Of particular interest, which we will include in our study, is the space of physical quantities, O.This will allow us to explore not only novel areas of the parameter space, but -perhaps more importantly -cover novel areas of the observables space, i.e. to explore novel phenomenological aspects of the model.To do so, we need to train a penalty p not on the values of the parameters, θ, but on its resulting physical quantities, {O i }, where the set {O i } can be of all or a subset of O i .
In our scans, we will include the novelty reward in both the parameter space and in the observable space, and in each case we will study penalties focused subsets of the parameters and the observables.

C. Further Implementation Details
We now discuss some implementation details of the ideas presented above.We implemented CMAES using DEAP -Distributed Evolutionary Algorithms in Python [45].The function C was slightly modified so as to force all values of C(O i ) to be nominally similar.
To achieve this, we have implemented in the code C(O i ) → log(C(O i ) + 1) which retains the same properties: positive semidefiniteness, continuous, and monotonically increasing away from the allowed interval.Furthermore, to prevent any constraint C(O i ) from dominating over any other, we rescale them at each generation to be bounded by [0, 1] using scikit-learn [46] MinMaxScaler before computing L and the final loss eq.(33).We used pyod [47] implementation of the HBOS [43], and set N bins = 100, observing a considerable computational overhead for higher values. 6he constraints that were checked in our main computational loop are listed in appendix A. For collider limits on novel scalars, we used HiggsTools version 1.1.3[35].As we perform the signal strength, µ ij , checks in our main computational loop, we only use the HiggsBounds functionality of HiggsTools, implemented using the HiggsBounds version 5 input files using the Python script provided by the HiggsTools authors, and used the HiggsBounds dataset version 1.4.
For each scan, we ran a total of 100 independent runs, each with a maximum of 2000 generations.We use the CMAES default parameters, which set the population size, λ, and the number of best candidates, µ, using a heuristic.The values for our problem were automatically set to λ = 12 and µ = 6, which means that each scan will at most evaluate 2000 × 12 × 100 = 2.4 × 10 6 points.CMAES has very few parameters to be defined by the user, only the initial mean of the multivariate normal and the overall scale of the covariance σ.For each run, we set the mean to a random point in the parameter space and initialise CMAES with σ = 1.
Furthermore, we follow the methodology in [27] where we define all operations related to CMAES and density estimation not over the parameter space, but over a hypercube [0, 1] d P , where d P is the number of parameters, that maps to the parameter space P.This allows us to treat all parameters in an equal nominal footing, avoiding any potential pathologies arising from having different parameters spanning many different orders of magnitude.
The final loss to be optimised to explore the parameter space is where p P ({θ j }) is the density penalty computed over the subset of parameters {θ j } effectively working as a novelty reward in P, p O ({O j (θ)}) is the density penalty computed over the subset of observables {O j } effectively working as a novelty reward in O, N C is the number of constraints.As discussed previously, we present different scans for different choices of {θ j } and {O j } to be included in the computation of p P and p O to promote focused scans.

VI. ANALYSIS AND RESULTS
We now present and analyse the results for multiple scans using the ideas presented in section V, and two random sampling scan strategies: purely random over the whole parameter space and 50% away from the alignment limit defined by eq. ( 27).We present two analysis, one where the HiggsBounds constraints using HiggsTools were not included in the loss function, and one where it has.The scan without HiggsTools runs faster in both computational overhead and CMAES convergence (to be discussed below), which allowed us to experiment with our methodology.The impact of HiggsTools on points obtained without including it in the loop is studied.We then perform a second analysis where we included HiggsTools in the loop to show how one can include multiple constraints from different sources and still be able to explore the whole parameter space of the model.
We start with scans that do not take into account HiggsTools results in the loss function.
All scans are performed over the whole 16-dimensional parameter space and all have the same 61 constraints.We then include HiggsTools in the optimisation loop by adding the respective contribution to the loss function.All scans and their details can be seen in table I.
As will be discussed in section VI D, HiggsTools reduces the number of successful runs by a factor of 2, and therefore for these runs we allow for 200 instead of 100 scans.We first study the implementation of CMAES with and without parameter space reward to show the enhanced exploration capabilities of CMAES when including a parameter density penalty in the loss function.In fig. 1 we show the scatter plot of the (sin(α 1 − β 1 ), sin(α 2 − β 2 )) plane for different runs.In particular, we exhibit the difficulty of random sampling finding valid points, with fig.1a showing only 23 valid points, of which only one passed the HiggsBounds constraints.These were obtained from a scan that sampled an estimated O(10 13 ). 7In fig.1b we show the points obtained by sampling within 50% of the alignment limit, where the allowed points are highly constrained with α 1 ≃ β 1 , an expected result due to the highly constraining bounds from collider measurements of the Standard Model-like Higgs boson decay channels.In the same graph, we can also observe the boundaries imposed by the alignment limit, as | sin(α 2 − β 2 )| ≃ 0.5.In the next plot, fig.1c, we show the result of the CMAES scan without further exploration.We observe a funnelling of the results into a single region α 1 , β 1 , α 2 , β 2 ≃ 0, clearly providing little coverage of the parameter space, although still providing far more valid points than the random sampler.This lack of exploration of CMAES was first observed in [27] and is easily understood from an algorithm point of view as CMAES does not have built-in mechanisms to escape a minimum (global or minimal). 8 7 We can only estimate the number of points as it would have been prohibitive to store all non-valid points.
Therefore, we measured how long the scan took to process a few thousand points, and kept a loose track of the random sampling run to produce the estimate. 8In [27] CMAES was endowed with a restart strategy, which mitigates this and allowed CMAES to draw a more complete picture of the valid regions of the parameter space.In this paper, we have not implemented this, as our focus is on developing a novelty reward driven exploration.A point surviving HiggsTools is represented in green, otherwise in red.
Despite the seemingly lacklustre result, the points in fig.1c obtained by CMAES were obtained in a quick run of only around O(10 3 ) attempts, providing around ten orders of magnitude sampling efficiency improvement over the random scan (see section VI D for a more detailed discussion on convergence metrics and performance).This allows us to change the sampling logic as to explore the parameter space once the sampler converges into a valid subregion of the parameter space.As explained in the preceding sections, this is achieved by including a parameter density penalty.In fig.1d we show the result for the CMAES runs when activating the novelty reward for all parameters, i.e. using a density penalty for all parameters.We can immediately see a much larger region of the parameter space scanned, especially beyond the alignment limit in the sin(α 1 − β 1 ) direction.We can also observe the first artefact of this methodology: we can see sequences of points, akin to the trail of a paintbrush, in this plane.These trails are in fact the path that CMAES has covered while exploring the parameter space away from previously found points.By introducing the parameter penalty over all parameter space, we were able to find novel points away from the alignment limit.However, because the penalty is computed over all the parameter space, CMAES has no incentive to explore interesting regions of the parameter space, as it can reduce the density penalty by spreading across parameters which have little impact on the constraints. 9To mitigate this, we focus the parameter density penalty on the four parameters described by these scatter plots, α 1 , β 1 , α 2 , β 2 .The resulting points can be seen in fig.1e, where we see that CMAES was able to spread its exploration even more in the (sin(α 1 − β 1 ), sin(α 2 − β 2 )) plane.More interestingly, the points that subsequently pass HiggsBounds, shown in green, cover a much larger region than those obtained using the sampling around the alignment limit, although the latter has arguably a better coverage over sin(α 2 − β 2 ) values more uniformly over different values of sin(α 1 − β 1 ).
The above scans were produced by performing a run without checking for the constraints coming from HiggsBounds provided by HiggsTools.The survival rate against HiggsTools of the points obtained using these scans is 3% − 5%, or, in other words, a factor of 1/20 or less.Furthermore, the execution time with HiggsTools increases by a factor of around 3.
To first approximation, checking for HiggsBounds in the loop can slow down the process of finding good points by an expected factor of 60 or more.On the other hand, as we can see in fig. 1 not using HiggsTools leads to too many invalid points, and depriving CMAES of this information will only prevent it from finding points that survive HiggsBounds.In fig. 2 we present the first scans with HiggsTools in the loop to check for HiggsBounds constraints.
The scans presented in fig.2a and fig.2b are direct analogous to those presented in fig.1d and fig.1e, respectively.In both cases we observe a far wider coverage of the parameter space than before, showcasing the importance of providing HiggsTools feedback to CMAES.
More importantly, both runs covered the space of valid points around the alignment limit in this plane, but went far beyond in the α 1 , β 1 subspace.We now turn to the masses of the charged scalars, which are constrained by direct searches.In fig. 3 we show the points obtained in the (m H + 1 , m H + 2 ) plane.The logic is similar to the previous discussion, with the difference that the last plot, fig.3e, shows the points from a scan where the penalty was focused on the charge scalar masses instead of the mixing angles.From the scan around the alignment limit, fig.3b, we can observe how HiggsBounds affects the valid region, especially for small values of scalar masses.Interestingly, in fig.3c we see that CMAES has provided more coverage over this cross section of the parameter space than before.This can be easily interpreted: CMAES works akin to a gradient descent algorithm, but with the performance enhanced by the approximation of local second derivative.This means that CMAES rolls down the loss function with momentum, following the quickest path to its minimum.This preference for a quick convergence path explains why different CMAES runs will provide similar values of the most constrained parameters, as it is through them that a path needs to be found so as to minimise the loss function.This eagerness to converge is a feature of CMAES, which is on the exploitive side of the exploration-exploitation trade off, as previously also discussed [27].In fig.3d and fig.3e we show the results of the scans with the parameter density penalty over all parameters and focused on the charged masses, respectively.We observe that both were able to cover more parameter space than the alignment limit random scan, but the scan focused on the charged scalar masses was able to cover more of the (m H + 1 , m H + 2 ) plane, especially in the m H + 1/2 ≳ 100 GeV limits.This result further shows that focussing the exploration reward on subsets of parameters can help uncover novel regions overlooked by traditional scans, although in this case most of the points with m H + 1/2 ≳ 100 GeV did not survive HiggsBounds.
Analogously to the discussion above on the mixing angles α 1 , β 1 , α 2 , β 2 , we now present the results with HiggsTools included in the loop to check HiggsBounds in fig. 4. We see that both with unfocused, fig.4a, and charged masses focused, fig.4b, novelty reward CMAES is able to cover a much larger parameter space region than the alignment limit sampling.Furthermore, the valid points found also span a larger region than those in fig. 3   The paths taken by CMAES while exploring the parameter space are very prominent in the scans just discussed.To better understand how CMAES is exploring, in fig. 5 we show the path traversed by a run projected onto the (m H + 1 , m H + 2 ) plane.This run has converged at generation number 129, with values (m H + 1 , m H + 2 ) ≃ (218.6,151.4) GeV.At generation 129, the overall scale of the covariant matrix, given by σ of the CMAES algorithm, is σ ≃ 0.002, a value much smaller than the initialised value of σ = 1.Once converged, the density penalty is then added to the loss function forcing CMAES to explore new values of the parameters, as can be observed in the left pane.As it explores, CMAES might be slowed down by the penalty; this leads the algorithm to increase σ to find new good points farther away.On the right pane we see this dynamical adaptation of σ by CMAES, with higher (lower) values of σ leading to more (less) localised samplings.This ability to adaptively increase σ when slowed down also provides CMAES with the capacity to escape local minima, and in our case it provides a way of forcing CMAES to move away from where it has been.∼ 140 GeV, which the scans with HiggsTools in the loop missed.We present a zoomed look of this region in fig.6, where we only show the points that have passed HiggsBounds constraints.This suggests that we have not completely mitigated the excessive eagerness of CMAES, which might lead us to miss multimodal solutions, i.e., disjoint valid regions of the parameter space.To better understand whether CMAES is being driven away from this region by its eagerness to converge, we performed a dedicated scan where we restricted the parameter space to m H + 1,2 ≤ 150 GeV, and with all other parameter bounds unchanged.We present the result in fig.7, where we notice that if restricted to that region, CMAES will explore it extensively.Furthermore, we notice that the points m H + 1/2 ∼ 140 GeV, which seem above to populate two disjoint regions, do not describe isolated islands of the valid parameter space.There are two important conclusions to draw from this.The first conclusion is that the empty regions of the scatter plot of valid points produced by CMAES do not equate to regions without valid points.This means that one has to be very careful when interpreting these seemingly empty regions without studying them in detail.The second conclusion is that when one focusses on studying these regions, one can find a completely different picture than assumed.In this case the previous results, both from alignment limit scan and from CMAES without HiggsTools in the loop, suggested that there are multiple disjoint regions of valid points in the (m ∼ 140 GeV, whereas a closer inspection teaches us that this is not the case.

B. Rewarding Exploration in the Observable Space
So far we have shown how we can improve the parameter space coverage by providing CMAES with a novelty reward in the parameter space by turning on a density penalty in observable values p({O i }).However, a more interesting avenue is to apply the novelty reward to the observable space, O,10 as this will allow us to assess whether there is new phenomenology obscured by traditional random sampling techniques.
In our first study we want to assess the impact on using a novelty reward in the observable versus the novelty reward in the parameter space studied above.In fig.8 we show the (µ ggF,γγ , µ ggF,ZZ ) scatter plots for different scans without HiggsTools in the loop.Similarly to the previous discussions on parameter space coverage, we see that CMAES without further exploration, fig.8c, provides a narrower coverage of the observable space than the alignment limit sampling strategy, adding to the intuition that CMAES alone is too eager to converge to be a reliable tool to draw a complete picture of the Physics.This changes considerably once we turn on the parameter space novelty award already studied, which also leads to a greater exploration of the observable space, as can be seen in fig.8d.This is easy to interpret, as forcing CMAES to explore the parameter space will always impact the value of the physical quantities of the model.We notice, however, that this is a byproduct of the parameter space exploration, as in this case CMAES does not have an explicit incentive to produce new observable values.In fig.8e we show the result of turning on the density penalty in the observable space, therefore explicitly forcing CMAES to find points with different phenomenology.The result is stunningly different from all the other scans, with CMAES finding points with a far more diverse value than any of the other scans considered so far.Of particular interest, we observe how CMAES can find points with µ ggF,γγ ≃ 1.2, which was completely obscured using alignment limit random sampling and painting a very different picture of what phenomenology the Z 3 3HDM model can have.Having shown how an observable space penalty can drive CMAES exploration into novel phenomenological realisations of the model, we now perform the scan with HiggsTools in the loop to endow CMAES with information of the HiggsBounds constraints (themselves added to the loss function and to the penalty).The resulting (µ ggF,γγ , µ ggF,ZZ ) scatter plot is shown in fig.9, where we observe how CMAES was able to find points across all allowed values (up to 2-σ with experimental measurement) for µ ggF,γγ with 0.89 ≲ µ ggF,ZZ ≲ 1.025, while completely rediscovering the possible values produced by the alignment limit sampling strategy.This result highlights the power and versatility of our methodology to find new phenomenological realisations of a model.µ Zγ , considering that, without additional states11 , the Higgs decay channels are considerably correlated, preventing any particular µ ij to be large while all others remain small.To study this, we performed a scan with a focused observable density penalty over (µ ggF,γγ , µ ggF,Zγ ), which we present in fig. 10 alongside the scatter plot obtained by the CMAES run with observable density computed over all constraints.Perhaps surprisingly, we see that the scan with the focused density penalty, fig.10b, has covered a smaller region than the one with the density penalty computed using all constraints, fig.10b.A possible interpretation for this is that the focused density is too constraining, preventing CMAES from finding other ways to populate this plane around other constraints.Conversely, the run with penalty over all constraints will be less demanding for CMAES to explore this subspace, as CMAES will find a way of reducing the penalty by spreading the possible constraint values elsewhere, eventually finding a new route to new values in the (µ ggF,γγ , µ ggF,Zγ ) plane.In other words, although when projected onto the (µ ggF,γγ , µ ggF,Zγ ) plane the valid region appears simply connected, the overall geometry and topology of the valid region of the parameter space are likely far more intricate with focused scans obfuscating these nuances.This interplay between a focused density, the availability of paths for CMAES to explore, and the topological and geometrical details of the valid region is an aspect of our methodology that will be further explored in future work.

C. Using Points as Seeds for New Runs
The scans performed so far have highlighted the versatility of our methodology in exploring parameter (and observable) spaces.However, the runs performed are independent of each other, i.e., while each run has its own parameter/observable density estimator, this is only trained using valid points found during that run alone.
Then, one can entertain the idea of reusing the information of previous scans to guide new runs in regions of interest.In this section, we explore this idea and provide an example of its implementation by choosing valid points from the previous scans as a seed to new runs.Recalling that CMAES can be initiated with an explicit mean, i.e. starting position, and an overall scale of the covariant matrix, σ, we can then use a valid point as the starting position of a new scan.In order to start exploring the vicinity of our starting position σ cannot be too large, and we found that setting it to σ = 0.01 guarantees that CMAES starts already at the minimum of the constraint loss function.
Seed points were identified by running HBOS on the entire collection of valid points (left pane of fig.11).For this concrete example, we evaluated the density only on the (m H + 1 , m H + 2 ) subspace and identified the 1% outliers (middle pane of fig.11), i.e., the points representing the least explored parts of the valid region of the parameter space.We notice some of the shortcomings of HBOS as the density estimator in this plot: given that HBOS fits a histogram to each dimension to compute the density, a point might be in a relatively sparse region but might not be picked as an outlier if one of its components is in a populated bin.
For example, we see that the outliers are not necessarily at the rim (convex haul) of the space, but in regions where there are few points both in m H + 1 and m H + 2 .This same shortcoming of HBOS is present in the scans with novelty reward, leaving room for improvement to be explored in future work.With the most outlying valid points identified, we ran 100 scans not only seeded by a point randomly chosen from the 1% outlier subset, but with the density penalty also making use of the outliers to guide the new scans away from the already explored regions.We did not use the whole sample of valid points to train the density estimator as it is comprised of over 4 million valid points, considerably slowing down the scan.In the right pane of fig.11 we show the resulting valid points found by the seeded scans, where we observe that CMAES was able to explore even further away from the previously chartered valid region.Clearly, one could now use the new points as new seeds in repeated iterations to explore even more this subsection of the parameter space, or any other section of it or of the observable space, in order to draw an even more global picture of the valid region.The caveat of only using chained scans is that one can only explore regions that are simply connected to the seed, a detail which must be kept in mind when employing this strategy.Finally, we notice that the scatter plot appears to have some vertical and horizontal regions with less points, this is an artefact of HBOS which draws a histogram with 100 × 100 bins in this subspace, impacting the density value along horizontal and vertical strips with width similar to the width of the bins.

D. Convergence Metrics
Having discussed how density penalties can be used to enhance the CMAES exploration of the parameter space, we now turn to another aspect of our methodology: the convergence speed.Recalling that CMAES operates by minimising the total loss function, L T from eq. ( 31), we show how its value decreases sharply after just a few generations in fig.12, where we also provide the values of random generations for comparison.More precisely, after just 100 generations, totalling around just 1200 points, CMAES has nearly converged to the valid region of the parameter space.Despite the suggestive previous plot, not all scans converge within the budget, which we set to 100 runs of 2000 generations for each case in table I without HiggsTools, and to 200 runs of 2000 generations for the cases with HiggsTools. 12We present these metrics in table II alongside statistics on the fraction of valid points that are within the alignment cases, AL-1 from eq. ( 26) and AL-2 from eq. ( 27).We see that, while most points are within AL-1, only a minority are within AL-2.More interestingly, we observe how the scan with novelty reward in the α 1 , β 1 , α 2 , β 2 subspace of the parameter space has produced the most points away from the alignment limit.This can be visually understood in fig. 1 and fig. 2 where it is clear that the novelty reward is guiding CMAES to values of α 1 , β 1 that are beyond the alignment limit bounds.This is a feature of the versatility of our methodology, as we can perform dedicated scans explicitly away from previously considered priors and regions of the parameter space. 12One could alternatively increase the budget for the HiggsTools by increasing the number of generations to 4000.Intuitively, more runs provide a more global picture, whereas longer runs allow for longer explorations of the valid region.The choice between allocating more budget to one over the other depends on the intended study.While the first CMAES rows correspond to scans without HiggsTools in the loop, the corresponding fractions of points in both alignment cases are computed using points that have passed HiggsBounds constraints after the scan.
As CMAES converges, it will start to find good points.This can be seen in fig.13 where we observe that after around 100 generations for the runs without HiggsTools in the loop, and after around 200 generations for the runs with HiggsTools in the loop, CMAES starts finding valid points.Table III: Sampling efficiencies of random sampling strategies.We note that for the completely random scan the numbers are estimated.
In fig.14  we see that the overall housekeeping overhead can be assessed in the random sampler rows as for these there is no overhead related to CMAES and to density estimation, and it is around 0.015 (0.019) seconds without (with) HiggsTools. 13The most important obser- 13 The larger housekeeping overhead associated with HiggsTools is due to the presence of more metrics to

VII. CONCLUSIONS
In this paper, we have developed a novel approach to explore the highly constrained multidimensional parameter space of the Z 3 3HDM, defined in section II and constrained keep track and larger intermediate files to save.
discussed in section III and appendix A, and go beyond alignment limit priors presented in section IV, by combining CMAES power of exploration with a Machine Learning estimator for point density.
It is important to note that, while the subject of study in this paper was the Z 3 3HDM parameter space, our approach is general and applicable to any Physics case, providing a solution to the difficulty of sampling good points in highly constrained multidimensional parameter spaces.
In section V we introduced our strategy, using CMAES, a powerful evolutionary strategy, in combination with HBOS, a fast ML model for density estimation.Our approach guarantees that the density-based novelty reward does not compete in the loss function with the constraints on the model and pushes CMAES to explore the parameter space once converged.Importantly, we showed how our methodology is versatile, as we can turn on the novelty reward in the parameter space or in the observable space, where the phenomenology is realised.Additionally, the novelty reward can be computed by estimating the density in only a subset of parameters and/or observables, allowing for quick focused scans on regions of interest.
In section VI we presented the results of multiple scans performed with our methodology, each with different combinations of parameters and/or observables on which the density penalty was computed.We showed how our approach can effortlessly go beyond the alignment limit sampling strategies, finding valid points in regions of the parameter space hitherto ignored by such sampling strategies.More precisely, using the novelty reward in the parameter space, both in whole or in a subset, in section VI A we have found that it is easy to go beyond the alignment limit in the (α 1 , β 1 ) plane.In the same analysis, we showed how our methodology also exposes regions of heavy scalar masses, even preferring it over the region excluded by HiggsBounds, which we explored in detail in section VI A 1 by restricting the parameter space to m H + 1,2 ≤ 150 GeV, finding a considerably different picture of that region of the parameter space that one would get from alignment limit sampling.
While we set ourselves to explore the parameter space with CMAES combined with a novelty reward, the Physics of the model resides its space of observables and physical quantities.In section VI B we have the results of scans where the density penalty was computed in the observable space instead of the parameter space.The results uncover novel possible phenomenological realisations of the 3HDM, an important contribution of this work that would not have been possible to achieve without our AI-based scan.In particular, we find that it is possible to accommodate Higgs decay signal strengths larger than one up to their current upper experimental bounds, a phenomenological signature not captured by alignment limit random sampling strategies.Given the versatility in exploring different observables, we set to study whether the Z 3 3HDM can explain the recent measurement of µ Zγ ≃ 2.2 by ATLAS and CMS [48], finding that µ Zγ ≲ 1.1 in the Z 3 3HDM, not a surprising result given that decay signal strengths are highly correlated in this model and one cannot get arbitrarily high values for one of them without spoiling the experimental measurements of the remainder (for other possibilities, see, for instance, the discussion in [49]).Finally, in section VI D we discuss the convergence metrics of our algorithm and compare it with the pure and alignment limit random sample strategies.Our methodology is orders of magnitude faster (both in number of tried points and wall time) than the random sampling strategies, providing a solution to the random sampling efficiency problem in highly constrained multidimensional parameter spaces.
Although the methodology presented in this work provides impressive speedup and efficiency improvement when compared to random sampling strategies, we have encountered some shortcomings and less appealing characteristics that we want to improve in future work.First, the methodology used in the analyses employs independent runs, each with its own density estimator from which the novelty reward is derived.An alternative approach is to share information across runs to ensure the novelty of the exploration.In section VI C we showed how such a strategy could be implemented, where we identified the valid region of the parameter space less populated (in the (m H + 1 , m H + 2 ) subspace) by our scans and then used some of the points in that region as a seed for new runs.The resulting new points were significantly different from the ones found previously, highlighting the potential for even further exploration by chaining runs together, a methodological detail that can be improved in the future.Second, in the same study, we encountered some artefacts arising from the binning nature of the HBOS, which can, in principle, be mitigated by using a different density estimator (or a different novelty detector).In our early exploration, we tried a variety of alternatives, all significantly slower than HBOS, making our methodology impractical.
Producing a better way of assigning the novelty reward could solve the binning problem and any manifestation of the curse of dimensionality produced by it.Third, we have observed that our methodology might not explore all possible regions as CMAES intuitively follows a path of fastest descent.This was particularly clear in section VI A 1 where we addressed the overlooked region of small charged scalar masses.By restricting the parameter space, we were able to populate that region easily, but the fact that it was not explored in the first place shows that we need to be careful when interpreting empty regions as regions without valid points.Lastly, we have observed that the geometrical and topological details of the valid region of the parameter space might impact the possible exploration paths of CMAES.This can have a profound effect on the results when there are disjoint, not simply connected, regions of the parameter space supporting good points.We leave to future work the development of a way to assess whether the scans are capable of capturing multimodal valid regions confidently.
Finally, our methodology opens up the possibility for a complete exploration of other N HDM (or any other BSM Physics) parameter spaces in light of the current highly constraining experimental results and theoretical conditions.As our work shows, this could lead to novel phenomenological realisations of these models and, ultimately, to the possibility of novel experimental signatures.We leave this phenomenological study for the future.

Figure 3 :
Figure 3: (m H + 1 , m H + 2 ) scatter plot for different runs without novelty reward.A point surviving HiggsTools is represented in green, otherwise in red.
that survived HiggsBounds constraints, highlighting the importance of including HiggsTools in the loop.

Figure 4 :
Figure 4: (m H + 1 , m H + 2 ) scatter plot for different runs with novelty reward and HiggsTools constraints included in the loss function.

Figure 5 :
Figure 5: Path of a CMAES scan with focused parameter density penalty on the (m H + 1 , m H + 2 ) plane.Left: colour representing the generation number.Right: colour representing σ, the overall scale of the covariant matrix of CMAES.
CMAES w/ constraint density penalty and HiggsTools in the Loop

Figure 11 :
Figure 11: (m H + 1 , m H + 2 ) scatter plot for the seeded run.Left: The whole collection of valid points obtained by the other scans.Middle: The 1% outliers classified by HBOS.Right: New points obtained by the seeded runs started at points randomly selected from the 1% outliers.

Figure 12 :
Figure 12: Total loss as function of generation.Only the first 500 generations are shown.The random sampler curves are over random generations of 12 poitns, the same population size as CMAES.The shaded regions represent 0.95 confidence intervals computed using a bootstrap of 100 runs.

Figure 13 : 4 ×
Figure 13: Number of valid points founds as function of generation.Only the first 500 generations are shown.The random sampler curves are over random generations of 12 points, the same population size as CMAES.The shaded regions represent 0.95 confidence intervals computed using a bootstrap of 100 runs.

Table I :
, β 1 , α 2 , β 2 focus 16 61 α 1 , β 1 , α 2 , β 2 None µ ggF,γγ and µ ggF,Zγ focus 16 61 None µ ggF,γγ , µ ggF,Zγ List of scans performed, where d P is the number of parameters, N C the number of constraints, p P the novelty reward over the parameter space, p O the novelty reward over the space of physical quantities O i .

Table II :
Convergence and coverage statistics of the different scans presented in table I.
we show the distribution of the generation where the first valid point was valid for the runs with and without HiggsTools in the loop.We see that checking for HiggsBounds postpones the discovery of the first good point by a factor of around 2 in terms of generation number.However, this is far better than with random sampling (both purely that CMAES without including HiggsBounds constraints in the loop tends to find points within O(100−500) seconds, i.e. in minutes, while when including HiggsBounds constraints, this increases do O(750 − 2250) seconds.The slowing down is easily understood: using HiggsTools slows down the evaluation of a point by a factor of around 3 (see below for Distribution of the elapsed time until the first valid point.The distributions shown are obtained from all runs with and without HiggsTools in the loop.To better understand the impact of the different components of our methodology in the total time, we present in table IV the times taken by different steps of the loop for CMAES and the random sampler, and with and without HiggsTools in the loop.In this table, generation time represents the time needed to perform all the steps of a generation including evaluation time, i.e. the time needed to compute all observables (including HiggsBounds, when applicable), train the density estimator (when applicable), and perform diverse housekeeping tasks such as save intermediate results, keep track of run metrics, etc.In this table

Table IV :
vation to take from this table is that the overall overheard of our methodology, including that associated with CMAES and density estimation, is at most around 10% of the total generation time for the CMAES runs without HiggsTools.Once we include HiggsTools in the loop, the overall generation time increases 3-fold but the overhead remains mostly the same, showing that our methodology provides even greater gains for problems with a slow evaluation time.Additionally, we see that our choice of HBOS for density estimator corresponds to a minor fraction of the overhead.However, one can notice that the standard deviation of the density estimator training is greater than the mean; this is because HBOS has a linear computational complexity with respect to the number of valid points, effectively becoming slower to train the more valid points we have found.Improvements to the density estimator are left for future work.Comparison of times taken by different parts of the loop for CMAES and random sampler with and without HiggsTools in the loop.The times refer to a generation of 12 points in every case.The values presented are the mean ± one standard deviation over the scans falling in each of the four categories.