Total: 1
We present the first systematic and comprehensive scan of two-flavon Froggatt-Nielsen (FN) models, employing artificial intelligence techniques to explore the high-dimensional, mixed discrete-continuous parameter space. Extending the standard single-flavon FN framework to a two-flavon setup in which separate flavon fields couple independently to the up- and down-type sectors, we demonstrate that the relative phase between their vacuum expectation values (vevs) provides a natural and generic source of CP violation absent in single-flavon models. To explore this enlarged model space, we cast the search for phenomenologically viable models as a multi-objective optimisation problem, formulating each experimental constraint as a separate objective, and employ the Non-dominated Sorting Genetic Algorithm III to simultaneously fit all 18 FN charges, 45 Wilson coefficients, and flavon parameters to both the quark and lepton sectors. Our approach requires no separate training phase and identifies phenomenologically viable models orders of magnitude faster than prior reinforcement learning methods. Imposing experimental constraints on CKM and PMNS mixing angles and CP phases, charged fermion masses, and neutrino squared-mass differences, we discover over $100\,000$ unique viable models with a remarkably low duplication rate, indicating that the space of valid two-flavon FN realisations has not been exhausted. Both Normal and Inverted neutrino mass squared orderings are realised, with the relative hierarchy between the flavon vevs producing qualitatively distinct predictions for the effective neutrinoless double beta decay mass $m_{ee}$. We further demonstrate the existence of minimal FN realisations with maximal flavon exponent as small as three, and of models reproducing charged fermion masses to within $6\%$ without any dedicated continuous parameter optimisation.