§ INTRODUCTION

The smallness of the observed neutrino masses (needed to explain neutrino oscillation experiments <cit.>)  is yet another consequence/success of the chiral nature of the Standard Model (SM) gauge structure. In the absence of, yet to be confirmed, right-handed neutrino at low energy, these neutrino masses cannot  emerge from marginal interactions among the SM degrees of freedom, like in the case for quarks and charged leptons, but require irrelevant interactions <cit.> possibly stemming from interactions with new physics degrees of freedom or new strong interactions (for a review, see Ref. <cit.>). A particularly compelling scenario is the type-I seesaw mechanism with right-handed sterile neutrinos endowed with a large SM gauge-invariant Majorana mass and coupled via Yukawa interactions to the left-handed active neutrinos <cit.>. The value of the scale of this Majorana mass, which also  measures the breaking of an accidental symmetry of the SM - the conservation of the total number of leptons, remains unknown and largely unconstrained.

In this paper, we will focus on the renormalizable theory of the SM extended with $n_N=3$ flavors of right-handed sterile neutrinos $N$, often referred to as or $\nu$MSM <cit.>. Depending on the mass we assign to these new particles, our discussion will capture both a type-I seesaw scenario and a scenario with light right-handed sterile neutrinos. Just like through the CKM matrix in the quark sector of the SM, adding more parameters to the theory by introducing a mass term opens the possibility to violate CP through physical phases that cannot be removed from the theory using field redefinitions. These new sources of CP violation (CPV) could help generating the matter-antimatter asymmetry via a leptogenesis mechanism <cit.>.

The purpose of this paper is the characterization of the flavorful parameters of the using the language of flavor invariants and in particular classifying the new sources of CPV that are introduced through the neutrino masses with Jarlskog-like flavor invariants <cit.>. In this analysis we will make use of the fact that – as physical quantities – all observables should in principle be expressible in terms of a finite number of flavor invariant objects. Such a set of invariants is referred to as the generating set or the basic set in the literature. This set enables us to express all observables as functions of them, which will allow us to learn about the structure of flavorful CP violation in the theory; for instance, at which order in the couplings certain phases can appear in observables. It is also straightforward to differentiate between the Dirac and Majorana nature of the mass term for the light sterile neutrinos in our formulation by taking the appropriate limits.

Furthermore, the set of invariants can be used to study the seesaw limit of the and make statements about the order in the effective theory at which certain UV parameters – in particular those which are CP-odd – can appear. For instance, we can show that while all parameters present in the effective theory of the appear suppressed with only one or two powers of the heavy Majorana mass, they can only appear in a flavor-invariant way with a suppression of 2 and 4 powers of the Majorana mass.

As a main result of our paper, we will construct the generating set of flavor invariants, which is the maximal set capturing all physical effects of the theory. For what concerns CPV effects, we will also endeavor to address the question of what is the minimal set that can capture the necessary and sufficient conditions of CP conservation (CPC) in all special spectra (for instance degenerate and vanishing masses) in the theory. In contrast to the quark sector of the SM, where all CPV effects are captured by a single CP-odd Jarlskog invariant, the introduces hundreds of CP-odd invariants in the generating set. This considerably complicates the analysis, making the identification of CPC conditions in the an ongoing challenge yet to be resolved.

In order to simplify the process of finding a generating set of flavor invariants, we will make use of tools from invariant theory that have been developed in recent years in the context of theoretical particle physics. Invariants have been used both in UV complete theories <cit.>
and effective field theories (EFTs) to characterize the parameters of the theory with respect to CP <cit.> and other symmetries of EFTs <cit.>. Some problems which are closely related to our analysis in this paper have been previously investigated in the literature. In particular, the with only two generations of charged leptons and neutrinos has been studied in Ref. <cit.> and the case of adding two generations of right-handed sterile neutrinos to the SM has been treated in Ref. <cit.>. There, the authors also use the flavor invariants to formulate the necessary and sufficient conditions for CPV in the model. The Hilbert series of the with three generations has previously been reported on in Refs. <cit.>. As we will see throughout this paper, generalizing the discussion from two to three generations is not an easy task. The more complex flavor group structure leads to a significantly more complicated algebraic ring structure of the invariants in the set, a so-called non-complete intersection ring. We will find that the number of generating invariants in the theory jumps from 38 in Ref. <cit.> for two generations of sterile neutrinos to 459 for three generations of sterile neutrinos, and there is an enormous amount of polynomial relations among them. As a consequence of the complicated algebraic structure of the non-complete intersection ring, there are non-trivial cancellations in the plethystic logarithm which is often used to count invariants given a set of building blocks and their transformation properties. Hence, one should take care when solely relying on these tools from invariant theory to build a complete basis of flavor invariants.

The paper is structured as follows. In Sec. <ref>, we will define more precisely the setup for neutrino masses that we consider in this paper, and introduce all tools from invariant theory we will need to build a generating set of flavor invariants. In Sec. <ref>, we will present our results for 3 generations of right-handed neutrinos including the Hilbert series, the generating set and the primary set. In Sec. <ref>, we will analyze the seesaw limit of the , where we discuss how the physical parameters of the full theory map to those in the effective theory. In Sec. <ref>, we will study the minimal set of CP-odd invariants used to capture the CPC conditions. We emphasize that this set should be analyzed without assumptions about the spectrum of the parameters. We show the results for the theory with 2 generations of fermions. The case of 3 generations turns out to be too hard to solve, but we nonetheless discuss an approach towards finding such a set. Finally, in Sec. <ref>, we draw some conclusions. In the Apps. (<ref>–<ref>), we provide additional information and supporting materials that complement the main body of the paper.

§ THE NUSM AND INVARIANT THEORY

§.§ SM extended with right-handed neutrinos

$SU(3)_L \times U(1)_L$     $SU(3)_e \times U(1)_e$     $SU(3)_N\times U(1)_N$     $U(1)_{L+e+N}$
$Y_e$    $\mathbf{3}_{+1}$    $\mathbf{\bar{3}}_{-1}$    $\mathbf{1}_0$     0

$Y_N$    $\mathbf{3}_{+1}$    $\mathbf{1}_0$    $\mathbf{\bar{3}}_{-1}$     0

$M_N$    $\mathbf{1}_0$    $\mathbf{1}_0$    $(\mathbf{\bar{3}}\otimes_s\mathbf{\bar{3}})_{-2}$     $-2$
$X_e=Y_e Y_e^{\dagger}$     $(\mathbf{3}\otimes\mathbf{\bar{3}})_{0}$     $\mathbf{1}_0$     $\mathbf{1}_0$     0

$X_N=Y_N Y_N^{\dagger}$     $(\mathbf{3}\otimes\mathbf{\bar{3}})_{0}$     $\mathbf{1}_0$     $\mathbf{1}_0$     0

$X_M=M_N M_N^{*}$     $\mathbf{1}_0$     $\mathbf{1}_0$     $(\mathbf{\bar{3}}\otimes\mathbf{3})_{0}$     0

The flavor transformation properties of the relevant Yukawa matrices and Majorana mass matrix treated as spurions. The subscripts of the $SU(3)$ representations denote the charge under the $U(1)$ part of the flavor symmetry group. Furthermore, $\otimes_s$ denotes the symmetric tensor product of the simple representations. The charges of all spurions under Abelian lepton number transformations are indicated in the last column. We also show the transformation properties of $X_e$, $X_N$ and $X_M$, which will be used in following sections.

Before we start with the analysis, we have to precisely define the theory we will work with. There are several ways to introduce a mass terms for neutrinos at low energies. In this paper, we will extend the SM particle spectrum by adding 3 generations of right-handed sterile neutrinos $N$, which in the literature usually goes under the name of .[In principle 2 generations of sterile neutrinos are enough to generate the observed neutrino masses at low energies <cit.>. As mentioned before this case has been treated in Ref. <cit.>.] The most general renormalizable Lagrangian that can be built with these fields and the symmetries of the SM gauge group is

\begin{eqnarray} \label{eq:L4}
\mathcal{L}_{\nu\text{SM}}&=&
\sum_{\psi}\bar{\psi}i \slashed{D}\psi
-\left[\frac{1}{2}\left(NCM_NN\right)+\bar{L}Y_NN\widetilde{H}+\bar{L}Y_e e H
+\bar{Q}Y_u u \widetilde{H}+\bar{Q}Y_d d H +\mbox{H.c.}\right]\nonumber\\
&-&\frac{1}{4}G^a_{\mu\nu}G^{a\mu\nu}-\frac{1}{4}W^I_{\mu\nu}W^{I\mu\nu} -\frac{1}{4}B_{\mu\nu}B^{\mu\nu}+(D_\mu H)^\dagger(D^\mu H)
-\lambda\left(H^\dagger H-\frac{v^2}{2}\right)^2,
\end{eqnarray}

where $\psi$ represents all fermion fields $\{Q,u,d,L,e,N\}$, $H$ is the Higgs doublet with the vacuum expectation value (vev) $v$, $\widetilde H_i=\epsilon_{ij}H^*_j$, and $G^a_{\mu\nu},~W^I_{\mu\nu},~B_{\mu\nu}$ are the field strengths of the SM gauge fields and $i, j$ are $SU(2)$ fundamental indices, $a=1,\dots,8$ is a color index, and $I=1,\dots,3$ is a weak isospin index. The gauge symmetries allow for the introduction of Majorana masses for the right-handed neutrinos $N$, which corresponds to a symmetric matrix $M_N$.[If lepton number is an exact symmetry of the classical Lagrangian (the symmetry is anomalous in the SM and the ), the Majorana mass term is forbidden, and the neutrino will have a Dirac mass. Then, the flavor structure of the lepton sector will have exactly the same form as the quark sector in the SM.] The $Y_{u,d,e,N}$ are the $3\times 3$ Yukawa couplings with the Higgs field. $C$ is the charge-conjugation matrix, and $D_\mu$ is the gauge covariant derivative.

One can easily check that the fermion kinetic term in $\cL_{\nu\text{SM}}$ is invariant under unitary $U(3)$ flavor transformations of the fermion fields. Assuming that this flavor symmetry is only softly broken by the Yukawa couplings and Majorana mass term, we promote all flavorful couplings to spurions under this symmetry, making the Lagrangian formally invariant. The corresponding transformation properties of the spurions under the non-Abelian part of the flavor group can be found in Tab. <ref>.

One immediate consequence of this assignment is that the presence of $M_N$ breaks lepton number, as the transformation properties of $M_N$ do not allow for non-trivial rephasings of $N$. Therefore, there will be additional physical Majorana-type phases in the spectrum of the theory. These phases can be removed in the SM because all interactions enjoy a symmetry under respectively rephasing the lepton and quark fields.

Later on in the paper it will prove useful to have an explicit parameterization for the flavorful matrices. One physical parameterization is given by[The detail about this physical parameterization can be found in Ref. <cit.>. However, they use left-handed Weyl fields to construct the theory, thus leading to a different convention compared to ours. Therefore, we provide the details focusing on the new convention in App. <ref>. In addition, we introduce new parameterizations that can be more conveniently used for the study of algebraic properties of the polynomial rings.]

\begin{equation}\label{eq:Parameterization}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\) \, , \quad Y_N = V \cdot \diag\(y_1,y_2,y_3\) \cdot W^{\dagger}, \quad M_N = \diag\(m_1,m_2,m_3\) \, ,
\end{equation}

\begin{equation}
\begin{split}
V =& U(\theta_{12},\theta_{13},\theta_{23},\delta) \cdot \diag\(1,e^{i \phi_1},e^{i \phi_2}\),\
W = \diag\(1,e^{i \phi_1^{\prime}},e^{i \phi_2^{\prime}}\) \cdot U(\theta_{12}^{\prime},\theta_{13}^{\prime},\theta_{23}^{\prime},\delta^{\prime}) \, ,
\end{split}
\end{equation}

and $U(\theta_{12},\theta_{13},\theta_{23},\delta)$ has been defined in Eq. (<ref>), which is a CKM-like matrix with a phase $\delta\in [0,2\pi)$ and three mixing angles $\theta_{ij}\in [0,\pi/2]$. $\phi_{1,2}\in [0,2\pi)$ and $\phi'_{1,2}\in [0,\pi)$ are additional phases. This parameterization correctly captures the 9 mass parameters, 6 mixing angles and 6 phases of the theory. A detailed counting from a symmetry perspective can be found in Tab. <ref>.

In the Dirac limit, the Majorana mass term $M_N$ is set to zero. Consequently, the mixing matrix $W$ and the phases $\phi_{1,2}$ become unphysical, as they can be absorbed by redefinition of the right-handed neutrinos. The flavor structure of the lepton sector then has a similar form as the quark sector. In this scenario, the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, which describes the mixing between different neutrino flavors, arises in analogy to the CKM matrix in the quark sector of the SM. In the Dirac case, only one Dirac CP phase is physical and characterizes the CP-violating effects in the lepton sector.

In the seesaw limit, the heavy Majorana fermions are integrated out, resulting in the emergence of light Majorana neutrino masses for the left-handed SM neutrinos at low energies. The mismatch between the diagonalization matrices in the charged lepton sector and the neutrino sector gives rise to a PMNS matrix. Due to the existence of the Majorana mass term, some rephasings are not allowed, and the two additional Majorana CP phases become physical.

The distinction between the Dirac and Majorana nature of neutrinos has important implications for the flavor structure and CP-violating phenomena in the lepton sector, which can be probed experimentally. While neutrino oscillation experiments, such as T2K and NOvA, have made significant progress in measuring the Dirac CP phase, the precise value of this phase is still not known with high accuracy <cit.>. Moreover, the existence of Majorana phases remains experimentally elusive. Neutrino oscillation experiments are primarily sensitive to the Dirac phase and do not probe the Majorana phases. The most promising avenue to probe the Majorana nature of neutrinos is through the search for neutrinoless double beta ($0\nu\beta\beta$) decay <cit.>, although it is important to note that even if observed, $0\nu\beta\beta$ decay would not directly measure the Majorana phases themselves.

Real Parameters     Imaginary Parameters     Total
$Y_e,Y_N,M_N$     $2\times 9+6$     $2\times 9+3$     45

$U_L,U_e,U_N$     $3\times 3$     $3\times 6-3$     24
Difference     15     6     21

Number of physical parameters in the generic vevs of the flavor spurions $Y_e,Y_N,M_N$ of the presented in Tab. <ref>. Note here, that $U_L,U_e,U_N \in SU(3)$ because the lepton numbers are broken by the presence of the vev of $M_N$ as discussed in the main text and cannot be used to remove imaginary parameters, hence the `$-3$' in the second column.

§.§ The Hilbert series and plethystic logarithm

In this section, we will briefly review useful tools from invariant theory, developed in Refs. <cit.> for operator bases and flavor invariants, that simplify the building process and characterization of the low-energy flavor invariants. We will mostly follow the notation introduced in Refs. <cit.> here. Our goal in this paper is to find a minimal set of flavor invariants that allows us to parameterize all observables in the theory in terms of those invariants.

The central object of this paper will be flavor invariants $\cI$ that are combinations of Lagrangian parameters of the theory invariant under the maximal possible flavor group of the renormalizable Lagrangian defined in Eq. (<ref>). In a first step, we want to find a set of invariants that allows us to express all remaining invariants in the theory as a polynomial of the invariants in the set. As mentioned above, this set of invariants is called the generating set or basic set $\cS_{gen}$, and always has finite cardinality for reductive groups <cit.>.

The generating set is a set of invariants, in which no invariant can be expressed as a polynomial of all other invariants in the generating set, i.e.,

\begin{equation}
\label{eq:genset}
\cI \neq P(\cS_{gen}\backslash \{\cI\})\,, \quad \forall \cI\in \cS_{gen}\,.
\end{equation}

If an invariant can be expressed as a polynomial of the generating invariants, it will be termed as an explicit relation in this paper.[Later on in this paper, we will also introduce the terms explicit dependency and explicit redundancy when such a relation is observed.]
But the invariants in the generating set may still be algebraically dependent, i.e., there could exist relations between them in the form of

\begin{equation}
P\(\cI_1,\dots,\cI_m\) = 0\,,
\end{equation}

which are called syzygies in the invariant literature. We want to stress here that these kind of polynomial relations do not have a linear term in any of the invariants. Otherwise, it would contradict the definition of the generating set provided in Eq. (<ref>).

Among the generating set of invariants there exists a set of invariants which are furthermore algebraically independent, they are the so-called primary invariants.[Note that this set is not unique and many different choices are viable as long as the set is algebraically independent. Once a set of primary and secondary invariants – which unlike the primary invariants still have algebraic relations among them – is chosen, all invariants in the theory can be expressed as a polynomial of them following the so-called Hironaka decomposition (see App. <ref> for more details).] The fact that they are algebraically independent implies that there exists no syzygy only comprising of invariants from the primary set. Another interesting result is that the number of physical parameters, i.e., the minimal number of parameters that are left after all transformations allowed by the symmetry group of the theory are used, is equal to the number of invariants in the primary set <cit.>.

A useful guide to construct those invariants is the so-called Hilbert series (HS)

\begin{equation} \label{eq:HSdef}
\cH(q) = \sum_{i=0}^{\infty} c_i q^i\,,
\end{equation}

which enumerates the number $c_i$ of all possible invariants that can be built from the given set of Lagrangian parameters labeled by $q$ at a given order $i$.

It can be shown that the HS can always be written as a fraction of two polynomials <cit.>

\begin{equation} \label{eq:NumDen}
\cH(q) = \frac{\cN(q)}{\cD(q)} \, ,
\end{equation}

where the numerator is of palindromic form, i.e., $\cN(q) = q^p \cN(1/q)$ with $p$ the highest power of $q$ in $\cN(q)$, and all terms in $\cN(q)$ come with a positive sign. The denominator is of the form $\cD(q) = \prod_{i=1}^{m} (1-q^{d_i})$, where the total number of factors $m$ counts the number of primary invariants corresponding to the physical parameters in the theory, while the exponents $d_i$ in each factor give the power of the spurion in the invariant. If the numerator is trivial, i.e., $\cN(q)=1$ and the complete generating set is given by the set of primary invariants, the ring is called a free ring.

Of course, most theories contain more than one coupling, and it can be convenient to count each coupling with its own spurion to simplify the identification of the invariants in the HS. For $n$ independent couplings in the theory that are used to build invariants, one defines the multi-graded HS

\begin{equation}
\cH(q_1,\dots,q_n) = \sum_{i_1=0}^{\infty}\dots\sum_{i_n=0}^{\infty} c_{i_1 \dots i_n} q_1^{i_1} \dots q_n^{i_n} \, ,
\end{equation}

where the coefficient $c_{i_1 \dots i_n}$ now count the number of invariants containing the spurions $\(q_1,\dots,q_n\)$ to the power $\(i_1,\dots,i_n\)$. We call these powers the degrees of the invariant while we call the sum of the degrees the order of the invariant.

Note that the multi-graded HS is no longer guaranteed to come in the form of Eq. (<ref>), which has a palindromic property in numerator with positive terms and the denominator counting the number of primary invariants. To still obtain this information, one can always take the single-graded limit of the HS,  $\cH(q_1,\dots,q_n) \to \cH(q,\dots,q)$, where all spurions in the theory are counted with the same parameter.

After describing the properties of the HS, we now introduce the mathematical methods for calculating it. One convenient way to calculate the HS for reductive Lie groups is the so-called Molien-Weyl formula, which for a single coupling transforming in the representation $R$ of the group $G$ is defined as

\begin{equation}
\cH(q) = \int d\mu_G \ \exp\(\sum_{k=1}^{\infty} \frac{q^k \chi_R\(z_1^k,\dots,z_d^k\)}{k} \) \equiv \int d\mu_G \PE\[\chi_R(z_1,\dots,z_d);q\] \, ,
\end{equation}

where $d\mu_G$ is the Haar measure of the group, $\chi_R(z_1,\dots,z_d)$ is the character of the representation $R$ of the group $G$ of rank $d$, and we have defined the plethystic exponential (PE) in the last step. There is a straightforward generalization of the Molien-Weyl formula for a multi-graded HS in a theory with several couplings transforming in different representations $R_i$

\begin{equation} \label{eq:MolienWeyl}
\cH(q_1,\dots,q_n) = \int d\mu_G \prod_{i=1}^n \PE\[\chi_{R_i}(z_1,\dots,z_d);q_i\] \,.
\end{equation}

To study some of the properties of the ungraded HS of a theory with several couplings, we can take the single-graded limit $q_i \to q$.

Another useful function is the so-called plethystic logarithm (PL) which is the inverse function of the PE that we just defined, i.e., $\PE^{-1}(f(x)) = \PL(f(x))$ and is defined as follows

\begin{equation} \label{eq:PL}
\PL\[f\(x_1,\dots,x_N\)\] = \sum_{n=1}^{\infty} \frac{\mu(n)}{n} \log\[f\(x_1^n,\dots,x_N^n\)\] \, ,
\end{equation}

where $\mu(n)$ is the so-called Möbius function.[The Moebius function is defined as

\begin{equation}
\mu(n) =
\begin{cases}
0 & \text{n has repeated prime factors} \\
1 & n=1 \\
(-1)^j & \text{$n$ is product of $j$ distinct prime numbers} \\
\end{cases} \, .
\end{equation}

In most cases the PL proves extremely helpful because we can simply read off the number of generating invariants and syzygies from the coefficients of the spurions at a given order in the spurions.[For instance, the PL of the quark sector of the SM can be easily calculated with the HS shown in Ref. <cit.>, and is given by $\PL(q)=2 q^2+3 q^4+4 q^6+q^8+q^{12}-q^{24}$. From the positive terms, it can be read off that there are two order-2 invariants, three order-4 invariants, four order-6 invariants, one order-8 invariant and one order-12 invariant, while the negative term shows there is a syzygy at order 24. Among these 11 generating invariants, 10 of them are algebraically independent, which map to the ten physical parameters of the quark sector, with an additional invariant capturing the sign of the CP phase.] In the case of a complete intersection ring,[A ring is classified as a complete intersection if the difference between the number of generating invariants and the number of syzygies is equal to the Krull dimension (the Krull dimension being the maximal number of algebraically-independent invariants). Otherwise, it is categorized as a non-complete intersection <cit.>.] where the PL is just a polynomial in the spurions, the positive terms can be identified with the generating invariants of the theory, and the negative terms correspond to the syzygies that exist among them <cit.>. However, for non-complete intersection rings, the PL becomes a non-terminating series, and it has been noted in the literature <cit.> that the leading positive terms, i.e., all positive terms up to the first term with a negative sign in the PL,[The terms in the PL is usually sorted according to their powers in the spurions. Refer to footnote <ref> for an example illustrating the PL of the quark sector.] can be identified with the generating invariants, and the leading negative terms, i.e., the first negative terms that appear after the leading positive terms correspond to syzygies.[All other terms in the non-terminating PL of a non-complete intersection ring after the leading negative terms have – to our knowledge – no meaning for the construction of a generating set beyond the fact that they appear in a special form of the HS, the so-called Euler form <cit.>. In this form the HS can be written as $\cH(q) = \prod_{n=1}^{\infty}(1-q^n)^{-b_n}$, where it can be shown that the $b_n$ are exactly the coefficients in the PL $\PL\[\cH\(q\)\] = \sum_{n=1}^{\infty} b_n q^n$.]

We will see later for the with 3 generations of sterile neutrinos, that this does not necessarily have to hold true, and the interpretation of the positive and negative terms in the PL have to be slightly changed for more complicated invariant rings. This has also been pointed out in the literature <cit.>. As we will see in the next section, this is also the main difference between the SM quark sector and the . While, the representations that the spurions of the SM quark sector transform in are still sufficiently simple to generate a complete intersection ring with a terminating PL, this is no longer true for the . Here, the representation the Majorana mass $M_N$ lives in complicates the ring structure significantly, leading to a non-complete intersection ring for both two and three generations of sterile neutrinos. This non-complete intersection ring can also simply arise when increasing the number of generations of fields in a theory. For instance, the complete intersection ring in the quark sector becomes non-complete when the quark has 4 generations <cit.>.

§ BUILDING AN INVARIANT BASIS FOR THE NUSM

§.§ Hilbert series of the nuSM

Before we start building invariants, we will first compute the HS and PL to set our expectations for the generating and primary set of invariants. We use the Molien–Weyl formula introduced in Eq. (<ref>) to calculate the HS for the with the spurions $Y_e, Y_N$ and $M_N$.[To get a dimensionless quantity we take $M_N$ to be divided by the only other mass scale in the problem, the Higgs vev $v$. Only then, one can compare invariants with a different number of insertions of $M_N$.] For that, we need the characters of the fundamental and anti-fundamental representations and the Haar measure of $U(3)$ which are given by <cit.>

\begin{equation}
\begin{split}
\chi_{U(3)}^{\mathbf{3}} & = z_1 + z_2 + z_3 \,, \\
\chi_{U(3)}^{\mathbf{\bar{3}}} & = z_1^{-1} + z_2^{-1} + z_3^{-1} \,,  \\
d\mu_{U(3)} & = \frac{1}{6!} \(\prod_{i=1}^3 \frac{dz_i}{2\pi i z_i} \) \(- \frac{\(z_2-z_1\)^2\(z_3-z_1\)^2\(z_3-z_2\)^2}{z_1^2 z_2^2 z_3^2} \) \,.
\end{split}
\end{equation}

From these, we can construct the characters for the representations of the flavorful Lagrangian parameters of the following Tab. <ref>. For instance, the character for $Y_N$ is given by

\begin{equation}
\chi_{Y_N} = \chi_{U(3)_L}^{\mathbf{3}}(z_1,z_2,z_3) \chi_{U(3)_N}^{\mathbf{\bar{3}}}(z_4,z_5,z_6) = \(z_1 + z_2 + z_3 \) \(z_4^{-1} + z_5^{-1} + z_6^{-1}\)  \,.
\end{equation}

The characters for all other spurions can be obtained in the same manner. Using the expression for the Molien–Weyl formula in Eq. (<ref>) with the same grading for all spurions, we can calculate the ungraded HS. The calculation involves the integral over the six variables $z_1,\dots,z_6$ over the contour $|z_i|=1$, which can be obtained by calculating the residues. The same calculation has been presented before, please refer to Refs. <cit.> for details. We find for the numerator of the ungraded HS

\begin{eqnarray}
\cN(q) &=& 1+q^{4}+5 q^{6}+9 q^{8}+22 q^{10}+61 q^{12}+126 q^{14}+273
q^{16}+552 q^{18}+1038 q^{20}\nonumber\\
&&+1880 q^{22}+3293 q^{24}+5441 q^{26}
q^{28}+13417 q^{30}+19867 q^{32}+28414 q^{34}+39351 q^{36}\nonumber\\
q^{38}+68220 q^{40}+85783 q^{42}+104588 q^{44}+123852 q^{46}
q^{48}+159328 q^{50}\nonumber\\
&&+173201 q^{52}+183138 q^{54}+188232 q^{56}+188232
q^{58}+183138 q^{60}+173201 q^{62}+159328 q^{64}\nonumber\\
+142559 q^{66}+123852
q^{68}+104588 q^{70}+85783 q^{72}+68220 q^{74}+52604 q^{76}+39351
&&+28414 q^{80}+19867 q^{82}
+13417 q^{84}+8712 q^{86}+5441
q^{88}+3293 q^{90}+1880 q^{92}+1038 q^{94}\nonumber\\
&&+552 q^{96}+273 q^{98}+126
q^{100}+61 q^{102}+22 q^{104}+9 q^{106}+5 q^{108}+q^{110}+q^{114},
\label{eq:num}
\end{eqnarray}

which has a palindromic form.[ Note that in the numerator Eq. (<ref>), a term $(1+q^2)$ can be factorized, and the same factor also appears in the denominator Eq. (<ref>).  If this factor were to be simplified, the HS would take a rational form with a numerator featuring some negative terms, in contradiction with the positivity requirement announced earlier. However, one can always multiply a factor $(1+q^k)^m$ in both numerator and denominator if there exists a factor $(1-q^k)^n$ (where $m\leq n$) in the denominator. This multiplication removes a factor of $(1-q^k)^m$ from the denominator while introducing a new factor of $(1-q^{2k})^m$. The total number of factors in the denominator does not change, and the numerator keeps its palindromic form with positive terms. This freedom indicates that there is ambiguity in determining the form of the HS if there is no further requirement of the HS. This is also reflected in the Hironaka decomposition, where the set of primary and secondary invariants are not uniquely determined. Please refer to the App. <ref> for details.]
The denominator is

\begin{eqnarray}
\cD(q) &=&   \left(1-q^2\right)^3 \left(1-q^4\right)^4 \left(1-q^6\right)^4
\left(1-q^8\right)^2 \left(1-q^{10}\right)^2 \left(1-q^{12}\right)^3
\left(1-q^{14}\right)^2 \left(1-q^{16}\right).
\label{eq:deno}
\end{eqnarray}

As expected, the powers of the factors in the denominator add up to 21, the number of physical parameters in the which is also the cardinality of the primary set. Our result of the ungraded HS is consistent with those found in Refs. <cit.>. We have furthermore calculated the multi-graded HS with different parameters $\{e,m,n\}$ counting the degrees of the couplings $\{Y_e,M_N,Y_N\}$, which we only show in App. <ref> due to its length. To obtain the results in this paper, we have developed our own  code that can efficiently calculate the Hilbert series. The code will shortly be published as a  package under the name  <cit.>.

Plugging the ungraded HS in Eq. (<ref>) to calculate the PL, we find furthermore

\begin{equation} \label{eq:HS_PL}
\begin{split}
\PL[\cH(q)] =& 3 q^2 + 5 q^4 + 9 q^6 + 10 q^8 + 19 q^{10} + 40 q^{12} + 66 q^{14} + 92 q^{16} + 70 q^{18} - 124 q^{20} \\
&- 703 q^{22} - 2039 q^{24} - 4391 q^{26} - 7472 q^{28} - 8522 q^{30} + 590 q^{32} + O(q^{34}) \, .
\end{split}
\end{equation}

We only show the PL up to order 32. For higher orders, both positive and negative terms will appear repeatedly in an infinite series, which implies that the has an algebraic structure of a non-complete intersection ring.

The usual interpretation of the PL suggests that the leading positive terms indicate a total of 314 generating invariants. However, our analysis reveals that this number is incorrect. Discrepancies begin to arise at order 16 in the PL, where the count of generating invariants exceeds 92. This is due to the non-complete intersection nature of the ring, resulting in non-trivial cancellations between the number of generating invariants and the number of syzygies. Concrete examples will be provided in the following sections to clarify this matter. It is worth noting that a similar cancellation was observed in a low-energy neutrino model in Ref. <cit.>, which also corresponds to a non-complete intersection ring.
There is no rigorous proof for this cancellation, however we can try to understand it as follows.

For a sufficiently simple invariant ring, the orders in the PL corresponding to the appearance of syzygies are well-separated from the orders corresponding to the appearance of the generating invariants. When the ring becomes more complicated, more invariants are needed to describe the full algebraic structure of the ring, hence there exist more generating invariants at higher orders. If syzygies containing the lower order generating invariants still appear at a similar order as in less complicated rings, there will be an overlap between the regions of positive and negative terms. This overlap will result in cancellations between the number of generating invariants and the number of syzygies. Therefore, one should be cautious when using the PL to count the number of generating invariants and the number of syzygies in a non-complete intersection ring. Observing a negative term in the PL does not necessarily imply the absence of generating invariants, but rather indicates the presence of more syzygies than generating invariants. Hence, the coefficient in the PL is the difference between the number of generating invariants and the number of syzygies.

Moreover, the coefficients in the ungraded PL at a specific order can be subject to cancellations from terms which have a different grading for the same total order in the multi-graded PL but cancel once the ungraded limit is taken. In this sense, we can not naively assume that the leading positive terms in Eq. (<ref>) can capture all generating invariants as we expected for the theory with a complete intersection ring. Hence, the multi-graded PL (see App. <ref>) is our main guide to check if we have found the correct number of generating invariants and syzygies at a given order in the spurions. It is also conjectured in Ref. <cit.> that the generating invariants are all captured by the terms prior to the pure negative order[In the multi-graded PL, we sort the terms according to the order (total degrees) of $[emn]$ (c.f. Eq. (<ref>)). In a given order, if all terms are negative, it is called a pure negative order.] in multi-graded PL. The pure negative terms occur at order 26 in Eq. (<ref>) in our theory, so the generating set should already be found within order 24. However, to test the conjecture, we also construct invariants up to order 26 to see that indeed no generating invariant can be found at this pure negative order. As we will see below, our analysis supports this conjecture.

§.§ Constructing the invariants

box/.style=rectangle, draw=black, thick, fill=white, text centered, minimum height=4em, align = center,
label/.style=text centered,align = center,

(box1) [box, minimum height=3em] All possible flavor invariants
516'101 invariants;

(box3) [box, below of=box1, yshift=-2cm]   Generating set of invariants $\cS_{\text{gen}}$
generates all invariants in the theory polynomially.
459 invariants;

(box3a) [box, minimum width=11.0em, minimum height=1.6em, below of=box3, yshift=-0.2cm, xshift=-2.27cm] 208 CP-even;
(box3b) [box, minimum width=11.em, minimum height=1.6em, below of=box3, yshift=-0.2cm, xshift=2.27cm] 251 CP-odd;

(box4a) [box, below of=box3a, xshift=-3.cm, yshift=-2.5cm] Primary set of invariants $\cS_{\text{prim}}$
captures physical parameters in the theory.
21 invariants;

(box4b) [box, below of=box3b, xshift=3.cm, yshift=-2.5cm] Minimal CPC set $\cS_{\text{min}}^{\text{CPC}}$
captures CPC conditions.
unknown number;

(box5) [box, below of=box3, xshift=0.cm, yshift=-7.5cm] Primary set of invariants in the seesaw limit

captures physical parameters in the EFT.

$\text{without power counting} \begin{cases}
\text{12 invariants}\sim\cS_{\{C_5\}}^{12} & (\text{includes only $C_5$ to any order)} \\
\text{21 invariants}\sim\cS_{\{C_5,C_6\}}^{21} & (\text{includes $C_5$ and $C_6$ to any order)}
\end{cases}$

$\text{with power counting} \begin{cases}
\text{3 invariants}\sim\tilde{\cS}_{4}^{3} & (\text{at mass dimension 4)} \\
\text{12 invariants}\sim\tilde{\cS}_{6}^{12} & (\text{at mass dimension 6)}\\
\text{21 invariants}\sim\tilde{\cS}_{8}^{21} & (\text{at mass dimension 8)}
\end{cases}$
[arrow] (box1) – (box3) node[midway,right] [label,xshift=0.2cm] Eliminating redundancies order by order
using relations of the form $I = P(I^\prime \in \cS_{\text{gen}})$;
[arrow] (box3a) – (box4a) node[midway,left,xshift=-0.7cm] [label] Denominator of HS;
[arrow] (box3b) – (box4b) node[midway,right,xshift=0.7cm] [label] Hilbert's Nullstellensatz;
[arrow] (box3) – (box5) node[midway,fill=white,yshift=1.4cm] Replace $M_N^*\to M_N^{-1}$ node[midway,fill=white,yshift=0.5cm] Identify $C_5,C_6, C_7$ substructures;

Flow graph of different invariant sets that appear in the analysis alongside the algorithms that are used to obtain one set from another. The number in the top box correspond to all single trace invariants up to a total order of 26 in all spurions. We find 208 CP-even and 251 CP-odd invariants that make up the generating set of the ring defined by the Lagrangian parameters and their transformation properties. A complete list of these invariants can be found in App. <ref>. The 21 algebraically independent invariants are selected from the CP-even generating set to form the primary set. These invariants accurately capture the 21 physical parameters of the theory. To determine the CPC conditions, the CP-odd invariants are selected from the CP-odd generating set. However, our program fails to find the minimal set due to the complexity of the theory. Detailed explanations on Hilbert's Nullstellensatz can be found in App. <ref>. In the seesaw limit, we replace $M_N^*\to M_N^{-1}$ for the generating invariants, and identify invariants with substructures of $C_5, C_6$ and $C_7$ as defined in Eq. (<ref>). The number of primary invariants with and without considering the total suppression of the invariants in the heavy Majorana mass is obtained by calculating the Jacobian rank of the identified invariants. A detailed analysis can be found in Sec. <ref>, where the primary sets $\cS$ and $\tilde{\cS}$ are shown explicitly.

Although the PL may provide some clues regarding the number of generating invariants, their specific form remains unknown. While it is possible to construct invariants for some simple models manually, in the case of complex models such as the , which involves hundreds of invariants in their generating sets, it becomes unfeasible to manually build them. Furthermore, as mentioned at the end of the last section, in a sufficiently complicated ring, the orders at which generating invariants and syzygies appear in the PL might overlap, hence leading to cancellations. To construct a generating set, one therefore cannot solely rely on the information provided by the PL.

Instead,  we will start by constructing all possible invariants up to a given order suggested by the conjecture on the PL, which is done by a graph-based method introduced in this section. The objective is to eliminate redundant invariants with explicit relations and create a generating set. This generating set will then be used to obtain a primary set. Furthermore, we aim to find the minimal set of invariants that can determine the CPC conditions based on the CP-odd generating invariants of the theory. The entire process is summarized in Fig. <ref>, and we will provide a detailed explanation of each step in the following sections.

[every circle node/.style=draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm,
every path/.style=-Stealth[length=2.5mm]]
[circle, label=[label distance=0.5mm]$e$] (e) at (0,0) 1;
[circle, label=[label distance=0.5mm]$L$] (L) at (2.5,0) 2;
[circle, label=[label distance=0.5mm]$N$] (N) at (5,0) 3;
[circle, label=[label distance=0.5mm]$N^*$] (Ns) at (7.5,0) 4;
[circle, label=[label distance=0.5mm]$L^*$] (Ls) at (10,0) 5;
[circle, label=[label distance=0.5mm]$e^*$] (es) at (12.5,0) 6;

(L) to[bend left=-45] node[above]$Y_e$ (e);
(e) to[bend left=-45] node[above]$Y_e^\dagger$ (L);
(L) to[bend left=-45] node[above]$Y_N$ (N);
(N) to[bend left=-45] node[above]$Y_N^\dagger$ (L);
(Ns) to[bend left=-45] node[above]$M_N$ (N);
(N) to[bend left=-45] node[above]$M_N^*$ (Ns);
(Ns) to[bend left=-45] node[above]$Y_N^T$ (Ls);
(Ls) to[bend left=-45] node[above]$Y_N^*$ (Ns);
(Ls) to[bend left=-45] node[above]$Y_e^*$ (es);
(es) to[bend left=-45] node[above]$Y_e^T$ (Ls);

[circle, label=[label distance=0.5mm]$d$] (d) at (0,-2.5) 7;
[circle, label=[label distance=0.5mm]$Q$] (Q) at (2.5,-2.5) 8;
[circle, label=[label distance=0.5mm]$u$] (u) at (5,-2.5) 9;
[circle, label=[label distance=0.5mm]$u^*$] (us) at (7.5,-2.5) 10;
[circle, label=[label distance=0.5mm]$Q^*$] (Qs) at (10,-2.5) 11;
[circle, label=[label distance=0.5mm]$d^*$] (ds) at (12.5,-2.5) 12;

(Q) to[bend left=-45] node[above]$Y_d$ (d);
(d) to[bend left=-45] node[above]$Y_d^\dagger$ (Q);
(u) to[bend left=-45] node[above]$Y_u^\dagger$ (Q);
(Q) to[bend left=-45] node[above]$Y_u$ (u);
(Qs) to[bend left=-45] node[above]$Y_d^*$ (ds);
(ds) to[bend left=-45] node[above]$Y_d^T$ (Qs);
(us) to[bend left=-45] node[above]$Y_u^T$ (Qs);
(Qs) to[bend left=-45] node[above]$Y_u^*$ (us);

The flavor invariant graph that can be used to construct all possible single trace flavor invariants in the . To any closed walk that follow the arrows, one can associate a single-trace invariant. Note that the graph for the SM would have a “holomorphic” structure, i.e.,  it has two separated branches involving separately only fields or only their conjugates. This changes in the , where the transformation properties of the Majorana mass $M_N$ connect the holomorphic and anti-holomorphic branches. More details can be found in the main text.

Flavor invariants from closed walks in a graph

The construction of flavor invariants in our paper heavily relies on the flavor invariant graph, which is inspired by Ref. <cit.>. We show the graph of the in Fig. <ref>, where the lepton sector and quark sector are presented in the top and bottom panels respectively. In a graph, the nodes represent the different fields and their conjugates, while the arrows are labeled with flavor matrices connecting various fields at the nodes. These graphs have the following two main advantages. First, they display the transformation rules for all flavor matrices. Taking the basic elements of the graph as an example,

[every circle node/.style=draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm,
every path/.style=-Stealth[length=2.5mm]]
[circle, label=[label distance=0.5mm]$F_i$] (A) at (0,0) $i$;
[circle, label=[label distance=0.5mm]$F_j$] (B) at (2.5,0) $j$;
(A) to node[above]$Y$ (B);

where the arrow labeled with a flavor matrix $Y$, starting from vertex $i$ to vertex $j$, which are labeled with the fields $F_i$ and $F_j$ respectively, indicating that the flavor matrix $Y$ should transform as $Y\to \neat{U_{F_i}^{} Y U_{F_j}^{\dagger}}$ under the flavor group. The numbers displayed in the vertex are simply for labeling purposes, and $F_i$ should be understood as the field at vertex $i$, not to be confused with the flavor index. For example, for the graph in Fig. <ref>, we can read off $Y_e\to U_LY_eU_e^{\dagger}$ and $Y_e^*\to \neat{U_{L^*}^{}Y_{e}^{*}U_{e^*}^{\dagger}=U_{L}^{*}Y_{e}^{*}U_{e}^{T}}$, where we have used the fact that $U_{f^*}=U_{f}^{*}$. Second, following the directional flow of the arrows, one can pass through vertices and edges, creating “paths”. In mathematical terminology, a path with repeated inclusion of vertices and edges is referred to as a “walk”. To every walk, one can associate an object with specific transformation under the flavor transformations. For instance, in the graph below, starting from vertex $i$ and following a sequence of arrows until reaching vertex $j$,

[every circle node/.style=draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm,
every path/.style=-Stealth[length=2.5mm]]
[circle, label=[label distance=0.5mm]$F_i$] (A) at (0,0) $i$;
[circle, label=[label distance=0.5mm]$F_k$] (P) at (2.5,0) $k$;
[circle, label=[label distance=0.5mm]$F_l$] (Q) at (6,0) $l$;
[circle, label=[label distance=0.5mm]$F_j$] (B) at (8.5,0) $j$;
(A) to node[above]$Y_1$ (P);
[dashed] (P) to node[above]$Y_2,\dots,Y_{n-1}$ (Q);
(Q) to node[above]$Y_n$ (B);

the product of flavor matrices $X\equiv Y_1\dots Y_n$ will transform as $X\to \neat{U_{F_i}^{} X U_{F_j}^{\dagger}}$. We can simply read off an example from Fig. <ref> that $\neat{Y_{e}^{\dagger} Y_{N}^{}\to U_{e}^{} Y_{e}^{\dagger} Y_{N}^{} U_{N}^{\dagger}}$. If there is a walk with $i=j$,
which means the walk is closed and ends at the starting vertex, we have $X\to \neat{U_{F_i}^{} X U_{F_i}^{\dagger}}$. This implies that we can construct a single trace invariant $\Tr(X)$ under the action of the group. As a result, flavor invariants are mapped to the closed walks in the graph, and one can identify a closed walk with each single-trace invariant in the theory.[One can also find other forms of the flavor invariants, but they can be converted to the single-trace forms. See Ref. <cit.> for the discussion.]

For instance, one can find the following correspondences between invariants and the closed walks in the quark and lepton sectors
(the red paths represent closed walks on the graphs denoted by the black arrows connecting the vertices)

\begin{equation}
\label{eq:walk}
\begin{split}
\begin{tikzpicture}[baseline={([yshift=-1ex]current bounding box.center)}, every circle node/.style={draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm},
every path/.style={-{Stealth[length=2.5mm]}}]
\node[circle, label={[label distance=0.5mm]95:$d$}] (d) at (0,-2.5) {{\footnotesize 7}};
\node[circle, label={[label distance=1mm]$Q$}] (Q) at (2.5,-2.5) {{\footnotesize 8}};
\node[circle, label={[label distance=0.5mm]80:$u$}] (u) at (5,-2.5) {{\footnotesize 9}};
\coordinate (dl) at ($(d) + (-0.3,0) $);
\coordinate (Qb) at ($(Q) + (0,-0.3) $);
\coordinate (ur) at ($(u) + (0.3,0) $);
\coordinate (Qu) at ($(Q) + (0,0.3) $);
\draw[rounded corners, red, style=](Qb) to[bend left=-45] (ur) to[bend left=-45] (Qu) to[bend left=-45] (dl) to[bend left=-45] (Qb);
\draw (Q) to[bend left=-45] node[above=3pt]{\textcolor{blue}{$Y_d$}} (d);
\draw (d) to[bend left=-45] node[above]{\textcolor{blue}{$Y_d^\dagger$}} (Q);
\draw (u) to[bend left=-45] node[above=3pt]{\textcolor{blue}{$Y_u^\dagger$}} (Q);
\draw (Q) to[bend left=-45] node[above]{\textcolor{blue}{$Y_u$}} (u);
\end{tikzpicture}&\sim 8\to 9\to 8\to 7\to 8\\
&\sim \Tr(\neat{Y_{u}^{}Y_{u}^{\dagger}Y_{d}^{}Y_{d}^{\dagger}})\,,\\
\begin{tikzpicture}[baseline={([yshift=-1ex]current bounding box.center)}, every circle node/.style={draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm},
every path/.style={-{Stealth[length=2.5mm]}}]
\node[circle, label={[label distance=1.0mm]110:$L$}] (L) at (2.5,0) {{\footnotesize 2}};
\node[circle, label={[label distance=1.5mm]90:$N$}] (N) at (5,0) {{\footnotesize 3}};
\node[circle, label={[label distance=0.5mm]70:$N^*$}] (Ns) at (7.5,0) {{\footnotesize 4}};
\coordinate (Ll) at ($(L) + (-0.3,0)$);
\coordinate (Nb) at ($(N) + (0,-0.2)$);
\coordinate (Nsr) at ($(Ns) + (0.3,0)$);
\coordinate (Nsr2) at ($(Ns) + (0.35,0)$);
\coordinate (Nu) at ($(N) + (-0.1,0.3)$);
\coordinate (Nl) at ($(N) + (-0.3,0)$);
\coordinate (Ll) at ($(L) + (-0.3,0)$);
\draw[rounded corners, red, style=](Ll) to[bend left=-50] (Nb) to[bend left=-45] (Nsr) to[bend left=-50] (Nl) to[bend left=-60] (Nsr2)to[bend left=-50] (Nu)to[bend left=-45] (Ll);
\draw (L) to[bend left=-45] node[above]{\textcolor{blue}{$Y_N$}} (N);
\draw (N) to[bend left=-45] node[above=2pt]{\textcolor{blue}{$Y_N^\dagger$}} (L);
\draw (Ns) to[bend left=-45] node[above=5pt]{\textcolor{blue}{$M_N$}} (N);
\draw (N) to[bend left=-45] node[above]{\textcolor{blue}{$M_N^*$}} (Ns);
\end{tikzpicture} &\sim 2\to 3\to 4\to 3\to 4\to 3\to 2\\
&\sim \Tr(\neat{Y_{N}^{}M_{N}^{*}M_{N}^{}M_{N}^{*}M_{N}^{}Y_{N}^{\dagger}})\,.
\end{split}
\end{equation}

Note that the walks are equivalently given as a chain of numbers corresponding to the vertices that are passed through. We only show here simple examples of walks in subsets of the graph. Obviously, one can consider more complicated walks involving more vertices and obtain new flavor invariants accordingly. Due to the cyclicity of the trace, the invariant associated to a closed walk is independent of the starting vertex, e.g. $ 8\to 9\to 8\to 7\to 8 = 7 \to 8 \to 9 \to 8 \to 7$ for the first example above. Based on this, we can rotate the numbers in the chain to arrange them in lexicographically smallest order. In addition, since the last vertex in a closed walk is always the same as the first vertex, we can omit the last number in the chain. To further shorten the notation, we can also remove the arrow, resulting in an integer representation of the walk notation. By following this approach, all single trace invariants are uniquely represented as integers. For instance, the two invariants above are represented by the two integers 7898 and 234343, respectively.

The inclusion of the Majorana neutrino mass $M_N$ creates an important difference between the lepton and the quark sectors.
The quark sector has two separate parts in the graph: the closed walks involve either the fields or their complex conjugate at the vertices, but never mix both. The walks are holomorphic or anti-holomorphic. The set of invariants built from the anti-holomorphic closed walks is equivalent of the one built from the holomorphic walks (see App. <ref> for details).
The introduction of the Majorana neutrino mass term $M_N$ complicates the scenario in the lepton sector. The two conjugate parts are connected through this new flavor matrix, making the invariant structure much more complex, as we will see below.

In the Dirac limit, $M_N\to 0$, the quark and the lepton sectors  obviously have the same flavor structure, and thus the flavor invariants have the same form. In App. <ref>, we show how to construct the flavor invariants in the quark sector based on the graph approach. The generating flavor invariants are shown explicitly, which can be easily mapped to the flavor invariants of the lepton sector in the Dirac limit.

For non-zero and finite $M_N$, we can systematically enumerate the closed walks corresponding to all single-trace flavor invariants up to arbitrarily high order. Based on the arguments of Ref. <cit.>, we are particularly interested in invariants up to order 26 which is the order of the first pure negative term in the multi-graded PL of the , and therefore it should be possible to obtain a set of generating invariants out of them. In App. <ref>, we present the brute-force algorithm we used to construct all the 516'001 single-trace invariants we obtain up to order 26. However, there is still a lot of redundancy among these invariants, and we can immediately reduce the set of invariants using some simple relations.

* Transpose redundancy

The flavor invariant associated with the walk $W_1\equiv ij\dots kl$ is always accompanied by another walk with primed vertices in reverse order $W_2\equiv l'k'\dots j'i'$, where $v'=7-v$ for the graph in the lepton sector, featuring a left and right mirror symmetry. The invariants generated from these two walks are the same, and they are connected through the trace identity $\Tr \big(A^T\big)=\Tr\big(A\big)$. For instance, $\Tr\Big(Y_eY_e^\dagger\Big)=\Tr\Big(Y_e^*Y_e^T\Big)$.

* Conjugate redundancy

The walk $W_1$ defined above is also associated with $W_3\equiv i'j'\dots k'l'$. The invariants generated by these two walks are conjugate to each other. Thus, both $\Tr(A)$ and $\Tr(A^*)$ will be generated in our construction. As the CP properties of $\Re\Tr(A)$ and $\Im\Tr(A)$ are more transparent than those of $\Tr(A)$ and $\Tr(A^*)$, we will trade $\Tr(A)$ and $\Tr(A^*)$, which are generated by the graphs, with $\Re\Tr(A)$ and $\Im\Tr(A)$, whenever a complex invariant is found.

* Cayley–Hamilton theorem

The Cayley–Hamilton theorem, along with its simple variations, enables us to eliminate invariants or establish relations among them. Note that the theorem is not fully utilized, as there are more complex variations that are not easy to employ. Therefore, we will implement a generic numerical algorithm to handle the remaining redundancies.

The details of these redundancies are discussed in App. <ref>. After eliminating these redundancies, we are left with a set of 8'666 invariants, which is still too large to form a set of generating invariants.

Construction of the generating set
After the set of invariants are pre-reduced with well-known trace relations, we search for explicit relations of invariants at a given order in terms of the other invariants at the previous lower orders. This allows us to identify the generating set of invariants which, by definition, does not have any explicit redundancy. To fully remove all explicit dependence and obtain a generating set, we introduce a numerical algorithm outlined in App. <ref>, which converts the problem of finding polynomial relations to a problem of solving finite system of linear equations. This algorithm has been used in different forms in Refs. <cit.> before. We have improved the algorithm for this work to avoid redundant syzygies[Redundant syzygies are those which have previously appeared in the algorithm at a lower order in the spurions and are multiplied by another syzygy, which has previously appeared, or some invariant of the generating set (or a sum of both), hence reappearing at a higher order. These kind of syzygies evidently do not carry any new information. This is also discussed in Ref. <cit.>, where the term “old relation” is used.] as we show in detail in App. <ref>. By using this method, we are able to generate all possible polynomial relations among the invariants at each degree, which includes both explicit relations and syzygies. Whenever an explicit relation is found, the corresponding invariant is removed from our set. When the algorithm has eliminated all explicit redundancies at a given order, it continues constructing relations at the next order including all the invariants that have survived the algorithm up to this point in its preliminary generating set, until the maximum order 26 is reached.

Running this algorithm up to order 26, the invariant list is reduced to a set without explicit dependencies. Our final set includes 459 invariants which form a generating set of the flavor invariants in the . Out of the 459 generating invariants 208 are CP-even and 251 are CP-odd. We want to stress again that all invariants in the theory can be captured by polynomials of these 459 flavor invariants. Hence, any observable in the theory is in principle expressible in terms of these invariants. In the following sections, we will further reduce the generating set to a primary set which captures all physical parameters in the theory. Additionally, we aim to reduce the CP-odd generating set to a minimal CPC set that captures all necessary and sufficient conditions for CPC in the  (c.f. Fig. <ref>). We provide the full generating set in App. <ref>, which is split into CP-even set and CP-odd set. From now on, we will use $\cI_i (\cJ_i)$ to represent the $i$th invariant in the CP-even (CP-odd) set.

Explicit relations and syzygies
To cross-check if the correct amount of generating invariants was found, we can use the information provided by the multi-graded PL. In order to match the numbers in the PL given by the number of generating invariants minus the number of syzygies at each order, we also have to find all syzygies at a given order. All of this can be achieved by following the process outlined above. Now, we will provide some examples for the polynomial relations we have found, which include both explicit relations and syzygies. Our program scans and checks all terms in the PL from lowest to highest order. Prior to order 12, our reduced invariant set accurately reproduces the terms in the PL. This means that no polynomial relation can be found, and all the invariants in our pre-reduced set are generating invariants. At degree $e^6n^6$ in the spurions, two invariants are found in our pre-reduced set, given by $\cI\equiv\Re\Tr(\neat{X_{N}^{2}X_{e}^{2}X_{N}^{}X_{e}^{}})$ and $\cJ_{10}\equiv \Im\Tr(\neat{X_{N}^{2}X_{e}^{2}X_{N}^{}X_{e}^{}})$, where $X_N$ and $X_e$ are defined in Tab. <ref>. From the multi-graded PL in Eq. (<ref>), a positive term $+e^6n^6$ is observed. According to the usual understanding of the PL, there should be only one generating invariant at this degree, which means one of the invariants $\cI$ and $\cJ_{10}$ is redundant. This is confirmed by our algorithm, we find that $\cI$ can be written as a polynomial of other lower degree CP-even invariants

\begin{equation}
\begin{split}
6 \cI =\  &\cI_3^3 \cI_1^3-\cI_3 \cI_6 \cI_1^3-3 \cI_3^2 \cI_7 \cI_1^2+3 \cI_3 \cI_{12} \cI_1^2-\cI_3^3 \cI_4 \cI_1+\cI_4 \cI_{11} \cI_1+\\
+&3 \cI_3^2 \cI_{13} \cI_1-3 \cI_3 \cI_{18} \cI_1+\cI_3 \cI_6 \cI_9-\cI_9 \cI_{11}+3 \cI_{12} \cI_{13}+3 \cI_7 \cI_{18}\,.
\end{split}
\end{equation}

With this explicit relation, the CP-even invariant $\cI$ becomes redundant, while no such explicit relation can be found for the CP-odd invariant $\cJ_{10}$; as a result the latter is included in our generating set. The same explicit relation was found in the discussion of the quark flavor invariants in Ref. <cit.>. The CP-odd invariant $\cJ_{10}$ in the lepton sector is analogous to the Jarlskog invariant $I_{6,6}^{(-)}$ in the quark sector as shown in Eq. (<ref>), and is the only CP-odd generating invariant in the that has no dependence on $M_N$. Note that $I_{6,6}^{(-)}=2i\,\Im\Tr(\neat{X_{u}^{2}X_{d}^{2}X_{u}^{}X_{d}^{}})$, while in our notation the additional factor “$2i$” is omitted.

Running our algorithm up to order 14, we get exactly the same number of generating invariants as indicated by the positive terms in the multi-graded PL. However, we run into a mismatch at degree $m^8n^8$. There are 10 invariants from the pre-reduced set, but we only find 9 explicit relations. As a result, we are left with one invariant $\cJ_{76}$ that can not be written as any polynomial of other invariants and by definition, it should be identified as a generating invariant. But looking at the multi-graded PL in Eq. (<ref>), one does not find a positive term at degree $m^8 n^8$ and would naively believe that there are no generating invariants at this degree. However, as we mentioned in Sec. <ref>, there can be non-trivial cancellations between the number of generating invariants and syzygies. Indeed, we also find another relation at degree $m^8 n^8$, which is given by

\begin{align*}
&~~~3 \cI_5^2 \cI_3^4-3 \cI_2 \cI_{10} \cI_3^4-8 \cI_2^3 \cI_8 \cI_3^3+12 \cI_2 \cI_5 \cI_8 \cI_3^3+8 \cI_8 \cI_{10} \cI_3^3+12 \cI_2^2 \cI_{16} \cI_3^3-24 \cI_5 \cI_{16} \cI_3^3+\\
&+30 \cI_2^2 \cI_8^2 \cI_3^2-18 \cI_5 \cI_8^2 \cI_3^2+36 \cI_{16}^2 \cI_3^2-6 \cI_5^2 \cI_6 \cI_3^2+6 \cI_2 \cI_6 \cI_{10} \cI_3^2+8 \cI_2^3 \cI_{14} \cI_3^2-12 \cI_2 \cI_5 \cI_{14} \cI_3^2+\\
&-8 \cI_{10} \cI_{14} \cI_3^2+10 \cI_2^3 \cI_{15} \cI_3^2-18 \cI_2 \cI_5 \cI_{15} \cI_3^2-4 \cI_{10} \cI_{15} \cI_3^2-60 \cI_2 \cI_8 \cI_{16} \cI_3^2-12 \cI_2^2 \cI_{20} \cI_3^2+\\
&+24 \cI_5 \cI_{20} \cI_3^2-24 \cI_2^2 \cI_{21} \cI_3^2+24 \cI_5 \cI_{21} \cI_3^2+24 \cI_2 \cI_{28} \cI_3^2-24 \cI_2 \cI_8^3 \cI_3+8 \cI_2^3 \cI_6 \cI_8 \cI_3+\numberthis\\
&-12 \cI_2 \cI_5 \cI_6 \cI_8 \cI_3-8 \cI_6 \cI_8 \cI_{10} \cI_3-24 \cI_2^2 \cI_8 \cI_{14} \cI_3+24 \cI_5 \cI_8 \cI_{14} \cI_3-48 \cI_2^2 \cI_8 \cI_{15} \cI_3+24 \cI_5 \cI_8 \cI_{15} \cI_3+\\
&+48 \cI_8^2 \cI_{16} \cI_3-12 \cI_2^2 \cI_6 \cI_{16} \cI_3+24 \cI_5 \cI_6 \cI_{16} \cI_3+24 \cI_2 \cI_{14} \cI_{16} \cI_3+48 \cI_2 \cI_{15} \cI_{16} \cI_3-16 \cI_2^3 \cI_{19} \cI_3+\\
&+24 \cI_2 \cI_5 \cI_{19} \cI_3+16 \cI_{10} \cI_{19} \cI_3+48 \cI_2 \cI_8 \cI_{20} \cI_3-48 \cI_{16} \cI_{20} \cI_3+48 \cI_2 \cI_8 \cI_{21} \cI_3-48 \cI_{16} \cI_{21} \cI_3+\\
&+48 \cI_2^2 \cI_{27} \cI_3-48 \cI_5 \cI_{27} \cI_3-48 \cI_8 \cI_{28} \cI_3-48 \cI_2 \cI_{39} \cI_3+3 \cI_5^2 \cI_6^2-6 \cI_2^2 \cI_6 \cI_8^2+6 \cI_5 \cI_6 \cI_8^2+6 \cI_2^2 \cI_{14}^2+\\
&-6 \cI_5 \cI_{14}^2+24 \cI_2^2 \cI_{15}^2-12 \cI_5 \cI_{15}^2-12 \cI_6 \cI_{16}^2+24 \cI_{20}^2+12 \cI_{21}^2-3 \cI_2 \cI_6^2 \cI_{10}+12 \cI_2 \cI_8^2 \cI_{14}+\\
&+24 \cI_2 \cI_8^2 \cI_{15}-10 \cI_2^3 \cI_6 \cI_{15}+18 \cI_2 \cI_5 \cI_6 \cI_{15}+4 \cI_6 \cI_{10} \cI_{15}+12 \cI_2 \cI_6 \cI_8 \cI_{16}-24 \cI_8 \cI_{14} \cI_{16}+\\
&-24 \cI_8 \cI_{15} \cI_{16}+24 \cI_2^2 \cI_8 \cI_{19}-24 \cI_5 \cI_8 \cI_{19}-24 \cI_2 \cI_{16} \cI_{19}-24 \cI_8^2 \cI_{20}+12 \cI_2^2 \cI_6 \cI_{20}-24 \cI_5 \cI_6 \cI_{20}+\\
&-48 \cI_2 \cI_{15} \cI_{20}-24 \cI_2 \cI_{14} \cI_{21}+6 \cI_2^3 \cI_{26}-6 \cI_2 \cI_5 \cI_{26}-12 \cI_{10} \cI_{26}-48 \cI_2 \cI_8 \cI_{27}+48 \cI_{16} \cI_{27}+\\
&+24 \cI_{14} \cI_{28}-24 \cI_2^2 \cI_{38}+24 \cI_5 \cI_{38}+48 \cI_8 \cI_{39}+24 \cI_2 \cI_{55}+12 \cJ_1^2=0\,,
\end{align*}

This polynomial relation indicates a redundancy of the CP-odd invariant squared $\cJ_{1}^2$ – not an explicit redundancy of an invariant – and therefore by definition it is a syzygy. Hence, combining the two findings we have

\begin{equation}
\PL\(e,m,n\) \supset + m^8 n^8 - m^8 n^8 \, ,
\end{equation}

explaining the non-trivial zero in the PL. At higher degrees in the spurions even stronger cancellations appear, that can even generate negative terms at orders where generating invariants may exist. Some examples are

\begin{equation}
\label{eq:cancellation_example}
\begin{split}
\PL\(e,m,n\) \supset &(21-1)\,e^4 m^4 n^8 + (3-1)\,e^2 m^4 n^{12}
+ (3-6)\,e^2 m^6 n^{10} +(2-6)\,e^2 m^8 n^{8}\, ,
\end{split}
\end{equation}

where we have used $(n_g-n_s)$ as a coefficient to indicate the $n_g$ generating invariants and $n_s$ syzygies at the corresponding degree.

Although there may be challenges when identifying generating invariants and syzygies, as long as the terms in PL are correctly interpreted, we are able to determine the correct number of the generating invariants and syzygies at each degree. Our algorithm accurately generates the terms in PL up to order 24 based on the updated description of the PL. However, at order 26, there are some mismatches due to “redundant syzygies”, which are not simply products of syzygies that appeared at a lower order in the spurions. We provide solutions to address these mismatches. For more details, please refer to App. <ref>. The order 26 is the first order that only has negative terms in the multi-graded PL. We have also confirmed that there are no further generating invariants at this order, as all the invariants are found to be a polynomial of lower order invariants.

Finally, our program successfully generates the correct number of generating invariants and syzygies at each degree of $[emn]$, reproducing the coefficients of the multi-graded PL shown in Eq. (<ref>). We have summarized the number of generating invariants and syzygies, as well as the CP-even and CP-odd counting at each degree in Tabs. <ref> and <ref>. In the new interpretation, the PL should be presented in a new form taking into account the non-trivial cancellations. The coefficient of each term should be replaced with $(n_g-n_s)$ as shown and described around Eq. (<ref>). According to the new multi-graded PL, one should be able to read off the correct number of generating invariants and syzygies at each degree easily. Due to its length, we will not show it here. Furthermore, we summarize the counting information at each ungraded order in Tab. <ref>. With this table, the ungraded PL shown in Eq. (<ref>) can also be revised to the new form

\begin{eqnarray}
\PL[\cH(q)] =&(3-0)\,q^{2}+(5-0)\,q^{4}+(9-0)\,q^{6}+(10-0)\,q^{8}+(19-0)\,q^{10}+(40-0)\,q^{12}\nonumber\\
&+(30-733)\,q^{22}+(5-2044)\,q^{24}+(0-4391)\,q^{26}+ O(q^{28})\,.\nonumber
\end{eqnarray}

One can find that up to order 26, the coefficients are exactly the same as those shown in Eq. (<ref>).

2|c|Order     $2$     $4$     $6$     $8$     $10$     $12$     $14$     $16$     $18$     $20$     $22$     $24$     $26$     Total
2*Generating     CP-even     $3$     $5$     $9$     $8$     $12$     $17$     $25$     $33$     $41$     $34$     $17$     $4$     $0$     $208$

CP-odd     $0$     $0$     $0$     $2$     $7$     $23$     $41$     $61$     $61$     $42$     $13$     $1$     $0$     $251$
2|c|Syz.     $0$     $0$     $0$     $0$     $0$     $0$     $0$     $2$     $32$     $200$     $733$     $2044$     $4391$     $7402$
2|c|PL$=$Generating$-$Syz.     $3$     $5$     $9$     $10$     $19$     $40$     $66$     $92$     $70$     $-124$     $-703$     $-2039$     $-4391$     $-6943$

The number of generating invariants and syzygies from order 2 to order 26, where the generating invariants are split into CP-even and CP-odd in the counting. The difference between the number of generating invariants and number of syzygies precisely aligns with the terms in the ungraded PL shown in Eq. (<ref>). In the last column, we list the total number of CP-even, CP-odd generating invariants, syzygies, and their difference. In the complete intersection ring, the difference between the number of generating invariants and number of syzygies should be the Krull dimension, which is 21 in our theory. The negative number shown here featuring a non-complete intersection ring.

At higher orders, the coefficients in the PL can no longer be connected to meaningful quantities like the number of generating invariants or the number of syzygies at a given degree <cit.>. For instance, at order 28, we can find two positive terms

\begin{equation}
\label{eq:PL_order28}
\PL\(e,m,n\) \supset +6 m^{14} n^{14} + 4 m^{16} n^{12} \, .
\end{equation}

However, all invariants constructed by brute force at these two degrees are redundant after applying the Cayley–Hamilton theorem. Therefore, there is no generating invariant, and these two positive terms must be misleading.[The only meaning we can give to the coefficients of these terms is due to the Euler form we mentioned in footnote <ref>.] The scanning program should terminate at the first pure negative order <cit.>, and we can be assured that we get a complete and minimal set of generating invariants. It is worth noting that identifying syzygies at each order is not necessary to obtain the generating set. The identification of explicit relations is sufficient for this purpose. Counting correct number of syzygies adds much complexity to our analysis, and it is only necessary to verify that we have found the correct number of generating invariants as indicated by the PL.

§.§ A primary set for the nuSM

To reduce the generating set of invariants to a primary set that captures all physical parameters of the theory, we can use the denominator of the HS as a guide.[In addition, the Hironaka decomposition (c.f. App. <ref>) documented in mathematical literature proves to be useful. The significance of the Hironaka decomposition lies in its ability to simplify the analysis of the invariants in the ring, particularly in scenarios where primary invariants are CP-even and secondary invariants are CP-odd. This enables a more focused investigation of the linear span of CP-odd secondary invariants, providing necessary and sufficient conditions for CP conservation <cit.>. Once a set of primary and secondary invariants is chosen, there will be a unique decomposition of any invariant in the theory in terms of the primary and secondary invariants in the form of Eq. (<ref>). This decomposition can be used to check if any invariant can be written in the form dictated by it, which indicates if the choice of primary and secondary invariants is a proper choice. However, according to Eq. (<ref>) and Eq. (<ref>), identifying all secondary invariants requires scanning the invariants up to order 114, which poses a significant challenge in our theory. Therefore, we will not attempt to present the invariants corresponding to the Hironaka decomposition.] The primary set is defined as a set of algebraically independent invariants. These invariants will capture all physical parameters in the theory. In order to reduce redundancies from the beginning, we choose the candidate sets with cardinalities equal to 21, which is the number of physical parameters in the case of the . One way to check if a candidate set is algebraically independent is to calculate the Jacobian with respect to all parameters in a given parameterization of the Lagrangian, for which we will use the parameterization from Eq. (<ref>). If the rank is equal to the number of physical parameters in the theory, a set of algebraically independent invariants is found. Following this procedure we find the following primary set of invariants

\begin{align*}
\label{eq:PrimSetnuSM}
\cI_{1}&=\,\Tr\left(X_e\right),~\cI_{2}=\,\Tr\left(X_M\right),~\cI_{3}=\,\Tr\left(X_N\right),~\cI_{5}=\,\Tr\left(X_M^2\right),~\cI_{6}=\,\Tr\left(X_N^2\right),\\
\cI_{7}&=\,\Tr\left(X_e X_N\right),~\cI_{8}=\,\Tr\left(Z_{MN}\right),~\cI_{9}=\,\Tr\left(X_e^3\right),~\cI_{12}=\,\Tr\left(X_eX_N^2\right),~\cI_{13}=\,\Tr\left(X_e^2X_N\right),\\
\cI_{15}&=\,\Tr\left(X_N Z_{MN}\right),~\cI_{23}=\,\Re\Tr\left( X_e X_N Z_{MN} \right),~\cI_{25}=\,\Tr\left(X_e^2 Z_{MN}\right),\\
\cI_{34}&=\,\Tr\left(\neat{X_{e}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),~\cI_{35}=\,\Tr\left(\neat{X_{e}^{} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),\numberthis\\
\cI_{47}&=\,\Tr\left(\neat{X_{e}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{N}^{*} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),~\cI_{50}=\,\Re\Tr\left(\neat{X_{e}^{} X_{N}^{} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),\\
\cI_{54}&=\,\Tr\left(\neat{X_{e}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),~\cI_{65}=\,\Re\Tr\left(\neat{X_{e}^{2} X_{N}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),\\
\cI_{79}&=\,\Tr\left(\neat{X_{e}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*2} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right),~\cI_{91}=\,\Re\Tr\left(\neat{X_{e}^{2} X_{N}^{2} Y_{N}^{} M_{N}^{*} M_{N}^{} M_{N}^{*} Y_{N}^{T} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}\right) \, ,
\end{align*}

where $\cI_i$ is the $i$th invariant in the CP-even set shown in App. <ref>, and we have defined $X_e= Y_e Y_e^{\dagger},X_N= \neat{Y_{N}^{} Y_{N}^{\dagger}},X_M=\neat{M_{N}^{} M_{N}^{*}}$ and $Z_{MN} = \neat{Y_{N}^{} M_{N}^{*} M_{N}^{} Y_{N}^{\dagger}}$. We want to stress here that the set of algebraically independent invariants is not unique. In particular, there are many sets that are compatible with the denominator of the ungraded HS and have a Jacobian rank of 21.

Furthermore, we have only chosen CP-even invariants from the generating set to form the primary set. However, it is also possible to include CP-odd invariants, provided they are algebraically independent. For instance, some CP-odd invariants are included in the primary set of the seesaw effective field theory in Ref. <cit.>. It might sound strange at first that CP-even invariants can be used to capture the physical phases in the theory. However, the primary invariants only capture the absolute value of the physical parameters. In particular, the sign of the phases in the theory are not described by the primary invariants. This is a well-known result in the quark sector of the SM <cit.>, where the Jarlskog invariant has to be added to the set of primary invariants to complete the generating set of invariants, enabling the expression of all observables in the theory in terms of flavor invariants.[As already mentioned in footnote <ref>, the square of the Jarlskog invariant, i.e., the square of the sign of the CKM phase, is in turn expressible in terms of CP-even invariants in agreement with the statement that the square of a sign is trivial.]

In addition, we have chosen the orders of the invariants in our primary set to follow those of the denominator of the HS in Eq. (<ref>). For instance, the first term $(1-q^2)^3$ of the denominator in Eq. (<ref>) indicates that there should be 3 invariants of order 2. However, following our discussion in footnote <ref>, one can change the numbers in the denominator of the HS by multiplying the numerator and denominator of the HS with the same factor $(1+q^k)^m$. In this case, more algebraically independent subsets of the generating set with cardinality 21 can function as a primary set.

§ THE SEESAW LIMIT

A phenomenologically interesting limit of the is the seesaw limit where $M_N$ is taken to be much larger than the electroweak scale $v$, allowing for an EFT description of the model. In this limit, the heavy neutrinos are integrated out, resulting in a matching between the low energy theory and the UV theory. The seesaw limit can be captured by the invariants in this paper by making the observation that $M_N^{-1}$, appearing in the EFT description, transforms in the same way as $M_N^*$. Hence, by replacing $M_N^* \to M_N^{-1}$ in all invariants, we can analyze this limit. After the replacement there are several structures appearing in the invariants which correspond to different orders in the EFT. In particular, we have <cit.>

\begin{equation} \label{eq:SeesawMatching}
\begin{split}
\frac{C_5}{\Lambda} & \sim \neat{Y_{N}^{} M_{N}^{-1} Y_{N}^{T}}\,, \\
\frac{C_6}{\Lambda^2} & \sim \neat{Y_{N}^{} (M_{N}^{*} M_{N}^{})^{-1} Y_{N}^{\dagger}} \,, \\
\frac{C_7}{\Lambda^3} & \sim \neat{Y_{N}^{} (M_{N}^{} M_{N}^{*} M_{N}^{})^{-1} Y_{N}^{T}} \,, \\
& \qquad\qquad\vdots \\
\end{split}
\end{equation}

where $C_5$, $C_6$ and $C_7$ are the Wilson coefficients of the Weinberg operator <cit.>, the operator $\cO_{Hl}^{(1)}$ in the Warsaw basis <cit.> and the operator $\cO_{lHD}^{(2)}$ in Ref. <cit.>, respectively.[Note, that at dimension 6 (and also at higher orders), there are other operators that can appear in the matching. $\cO_{Hl}^{(1)}$ at dimension-6 and $\cO_{lHD}^{(2)}$ at dimension-7 are only examples of an operator that receives a contribution from the matching of the to the SMEFT <cit.>. However, our choice of Wilson coefficients ensures the lowest overall suppression of the invariants with the Majorana mass.]

As before, it is insightful to construct the invariant graphs to understand more easily which kind of structures can appear in the flavor invariants. The relevant graph for this limit is show in Fig. <ref>. Here, we show the graph for two types of the Wilson coefficients $C_{LL^T}$ and $C_{LL^{\dagger}}$ as an example, where $C_{LL^T}\to U_LC_{LL^T}U_L^T$ and $C_{LL^\dagger}\to U_LC_{LL^\dagger}U_L^\dagger$ respectively. For instance, $C_{LL^T}\sim C_5, C_7$ and $C_{LL^{\dagger}}\sim C_6$. This shows that the graph approach serves as a general method for constructing flavor invariants and can be employed effectively in both UV and EFT studies. Furthermore, the graph in the EFT allows for a clearer understanding of the structure of the invariants. For instance, when constructing invariants for the EFT with the Wilson coefficient $C_5$, a closed walk in the graph indicates that $C_5^*$ must also be included. Following this conclusion, it becomes clear that all constructed invariants must adhere to an even power counting.

[every circle node/.style=draw, fill=blue!10, inner sep=0pt, minimum size=3.5mm,  every path/.style=-Stealth[length=2.5mm]]
[circle, label=[label distance=0.5mm]$e$] (e) at (0,0) 1;
[circle, label=[label distance=0.5mm]$L$] (L) at (2.5,0) 2;
[circle, label=[label distance=0.5mm]$L^*$] (Ls) at (5,0) 3;
[circle, label=[label distance=0.5mm]$e^*$] (es) at (7.5,0) 4;

(L) to[bend left=-45] node[above]$Y_e$ (e);
(e) to[bend left=-45] node[above]$Y_e^\dagger$ (L);
[draw=green!70!black, loop above] (L) to [out=55, in=125, distance=2cm] node[above] $C_6$ (L);
[draw=green!70!black, loop above] (Ls) to [out=55, in=125, distance=2cm] node[above] $C_6^*$ (Ls);
[red] (L) to[bend left=-45] node[above]$C_5,C_7$ (Ls);
[red] (Ls) to[bend left=-45] node[above]$C_5^*,C_7^*$ (L);
(Ls) to[bend left=-45] node[above]$Y_e^*$ (es);
(es) to[bend left=-45] node[above]$Y_e^T$ (Ls);

The flavor invariant graph that can be used to construct all possible single trace flavor invariants in the seesaw limit of the , where $C_5$, $C_6$ and $C_7$ are the Wilson coefficients of the Weinberg operator <cit.>, the operator $\cO_{Hl}^{(1)}$ in the Warsaw basis <cit.> and the operator $\cO_{lHD}^{(2)}$ in Ref. <cit.>, respectively, and they transform as $C_{5,7}\to \neat{U_{L}^{} C_{5,7}^{} U_{L}^{T}}$ and $C_6\to \neat{U_{L}^{} C_{6}^{} U_{L}^{\dagger}}$. We show in the main text that it is sufficient to consider $C_5$, $C_6$ and $C_7$ to capture all parameters of the full theory.

We can make several interesting observations using our generating invariants. Keeping the invariants with only the structures of $C_5$ and $Y_e$ and calculating the Jacobian of the invariants with respect to all parameters in a chosen parameterization, we find that its rank is 12. This corresponds to the six masses, three mixing angles and three phases that appear in the low-energy theory of three flavors of neutrinos with a Majorana mass term $\cL \supset - 1/2 \, \bar{\nu}_L m_{\nu} \nu_L^c + \text{H.c.}$. There are a total of 15 invariants with only $C_5$ and $Y_e$ structure, we show the 12 algebraically independent invariants as follows:

\begin{equation}
\cS_{\{C_5\}}^{12} \equiv \{\cI'_{1},\, \cI'_{4},\, \cI'_{9},\, \cI'_{14},\, \cI'_{22},\, \cI'_{34},\, \cI'_{35},\, \cI'_{54},\, \cI'_{79},\, \cJ'_{101},\, \cJ'_{168},\, \cJ'_{221}\}\,,
\end{equation}

where $\cI'_i(\cJ'_i)$ can be obtained from our generating invariants $\cI_i(\cJ_i)$ (refer to App. <ref>) by the replacement $M_N^*\to M_N^{-1}$. We introduce the notation $\cS_{c}^{r}$ to indicate that the EFT only involves the Wilson coefficients in the list $c$, and the corresponding theory has a rank-$r$ primary invariant set.

By also keeping the invariants with the structures corresponding to $C_6$ in Eq. (<ref>), we find that the Jacobian of the invariants has rank 21, corresponding to the number of physical parameters in the . We choose the primary set to be

\begin{equation}
\cS_{\{C_5,C_6\}}^{21} \equiv \cS_{\{C_5\}}^{12} \cup \{\cI'_{8},\, \cI'_{17},\, \cI'_{25},\, \cJ'_{20},\, \cJ'_{58},\, \cJ'_{62},\, \cJ'_{124},\, \cJ'_{125},\, \cJ'_{189}\}\,.
\end{equation}

This has previously been noted in Ref. <cit.>, where this analysis is performed in the effective theory of the seesaw model. They claim that the physical information about the full theory can be completely covered by the EFT with only two operators corresponding to the Wilson coefficient $C_5$ and $C_6$. Here, we come to the same conclusion in the full theory of the .

Note that, up to now we have not considered the power counting of the effective theory, i.e., we also allowed for invariants with several insertions of the $C_5$ structure to obtain the ranks of 12 and 21. On the other hand, working in the effective theory by taking the seesaw limit, one should work consistently to a certain order in the power counting of the theory, which in this case is defined by the characteristic scale of the process divided by the Majorana mass. By counting the number of insertions of $M_N$, we find that one cannot reach rank 12 and 21 at dimension 5 and 6, respectively. Instead, without any insertions of $M_N$, one can reach rank three, corresponding to the masses of the charged leptons. The primary set is simply formed by three flavor invariants that only involve $Y_e$:

\begin{equation}
\tilde{\cS}_{4}^{3} \equiv \{\cI'_{1},\, \cI'_{4},\, \cI'_{9}\}\,,
\end{equation}

where $\tilde{\cS}_{d}^{r}$ is the rank-$r$ primary set of invariants at mass dimension $d$. At mass dimension 6 (with 2 insertions of $M_N$), rank 12 can be reached, and the corresponding primary set is given by:

\begin{equation}
\tilde{\cS}_{6}^{12} \equiv \tilde{\cS}_{4}^{3} \cup \{\cI'_{8},\, \cI'_{14},\, \cI'_{17},\, \cI'_{22},\, \cI'_{25},\, \cI'_{34},\, \cI'_{35},\, \cI'_{54},\, \cI'_{79}\}\,.
\end{equation}

At dimension 8 (with 4 insertions of $M_N$), one can reach rank 21, which correspond to the primary set:

\begin{equation}
\tilde{\cS}_{8}^{21} \equiv \tilde{\cS}_{6}^{12} \cup \{\cI'_{16},\, \cI'_{21},\, \cI'_{24},\, \cI'_{32},\, \cI'_{37},\, \cI'_{51},\, \cI'_{53},\, \cI'_{78},\, \cI'_{112}\}\,.
\end{equation}

This set also includes invariants with the structure of $C_7$. Hence, even though all the information about the is in principle contained in the Wilson coefficients $C_5$ and $C_6$, this information can only be accessed at higher mass dimensions. Considering the power counting of the theory, full rank is first reached at dimension 8 by also considering the structure of $C_7$ in the invariants instead of going to higher mass dimensions beyond dimension 8 by only including the structures of the Wilson coefficients $C_5$ and $C_6$. We have summarized the results in Fig. <ref>.

§ CONDITIONS FOR CP CONSERVATION

The invariant structures of the have been presented in various papers for different numbers of lepton flavors. In Ref. <cit.>, the generating set for the with only two generations of right-handed neutrinos and charged leptons ($n_N=n_f=2$) has been shown, but the discussion of CPC conditions has not been expanded on. In Ref. <cit.>, the minimal seesaw model with two generations of right-handed neutrinos and three generations of charged leptons ($n_N=2, n_f=3$) has been discussed. However, the discussion of CPC conditions is provided only with some assumptions about the parameter spectrum. For the three-generation case ($n_N=n_f=3$), the discussion of CPC conditions is still lacking, as the generating set has not been constructed, although the Hilbert series has been calculated in Refs. <cit.>. In Ref. <cit.>, without using the Hilbert series and explicit construction of the generating set, six CP-odd invariants are provided to characterize CPV effects of the three generation case. Also there, assumptions are made regarding the spectrum of the theory.

The spirit of using CP-odd flavor invariants to determine CPC conditions is that they can capture all possible CPC conditions of the theory, which not only include the vanishing phases, but also conditions that lead to unphysical phases, such as degenerate masses and texture zeros. This paradigm is well-established in the quark sector of the SM, where the Jarlskog invariant captures all such conditions leading to the CPC of the SM. In this paper, our objective is to identify a set of Jarlskog-like invariants capturing all possible CPC conditions of the .

It should be noted, that without the guidance of the invariants, it is not easy to derive all CPC conditions of a theory, as some special conditions may appear that are beyond the conditions as we mentioned above. We show the complete list of CPC conditions for the with $n_N=n_f=2$ in App. <ref>, where in one particularly interesting case a scenario with pseudo-real couplings[By pseudo-real couplings, we denote a set of couplings which have irremovable phases but still conserve CP. This is the case in models with discrete symmetries, when the effect of a CP transformation on the Lagrangian parameters can be undone by a field redefinition, even if there exists no basis where all couplings can be made real. This has been previously noted in the context of Three Higgs Doublet Models <cit.> and toy models with more involved discrete symmetry groups <cit.>. We show an explicit example for the with two generations in App. <ref>.] is found. As a result, studying CPC conditions through the generating CP-odd invariants provides a more reliable approach. In the following, we will define the minimal invariant set that captures the CPC conditions, and propose useful methods to find it. However, due to the complexity of the invariant structure of the , finding the proposed minimal set in general can be challenging. The main goal of this section is to establish a framework for studying the CPC conditions, with the final solutions to these sets left for future work.

The CP-odd invariants in the generating set $\cS_{\text{gen}}$ act as generators of all CPV observables. When these invariants vanish, it establishes the necessary and sufficient conditions for CP conservation. However, this does not mean that the CP-odd generating set is the minimal set characterizing this property. The CP-odd generating set is required to generate any value of the CP-odd invariants in a given parameterization, while the set that can determine if CPV exists only captures their roots. In principle, the latter set should be a subset of the CP-odd generating set, and from now on, we will call it the CPC set. If the invariants in the CPC set is reduced to the minimal number, the CPC set becomes the minimal CPC set $\cS_{\text{min}}^{\text{CPC}}$.

By definition, the minimal CPC set is sufficient enough to generate all roots corresponding to the CPC conditions, all other CP-odd invariants should vanish automatically on these roots. It is possible to find a minimal CPC set, such that

\begin{equation}\label{eq:CPVsetDef}
\cJ_{\text{min}} = 0, \, \forall \cJ_{\text{min}} \in \cS_{\text{min}}^{\text{CPC}} \Longrightarrow \cJ = 0, \, \forall \cJ \in \cS_{\text{gen}} \backslash \cS_{\text{min}}^{\text{CPC}}.
\end{equation}

The straightforward method one can try is to solve the common zeros of polynomials in a candidate minimal CPC set,[We have done this for the generating set of the ($n_N=n_f=2$) with the package  <cit.>, for which the complete set of CPC conditions are listed in App. <ref>. However, the common zeros of the subset of the generating invariants are still difficult to solve, thus we can not determine the minimal CPC set.] and subsequently apply the solutions to the other CP-odd invariants to check whether they will vanish. However, this is not practical for complicated polynomials. Without directly solving the polynomial equations, one can also use the syzygies to determine whether other invariants are automatically zero given that all invariants in a minimal set are set to zero. This approach was e.g. followed in Ref. <cit.>, where the author found some syzygies that can determine the minimal CPC set in the 2HDM. In this section, we will show that there is a specific form of the syzygy that can help to determine whether an invariant will vanish.

The general problem concerning the common zeros of polynomials is closely connected to Hilbert's Nullstellensatz <cit.>, a theorem that establishes a fundamental relationship between geometry and algebra. We have presented this theorem and relevant mathematical terms in App. <ref>. In this section, we will employ Hilbert's Nullstellensatz to reframe the problem of identifying the minimal CPC set. In the invariant ring, or more generally in the polynomial ring $R:=\mathbb{Q}[x_1,\dots,x_n]$, where $x_{1,\dots,n}$ are the parameters in the theory. The CP-odd generating invariants are polynomials in this ring. Let $\cJ_s,\cJ_1,\dots,\cJ_m\in R$, Hilbert's Nullstellensatz says that if an invariant $\cJ_s$ vanishes on all the common zeros of the $\cJ_{1,\dots,m}$, then there exist some integer $t$, such that $\cJ_s^t$ is a subset of the ideal $I:=(\cJ_1,\dots,\cJ_m)$, i.e.,

\begin{equation}
\label{eq:jpower}
\cJ_s^t=f_1 \cJ_1+\dots+f_m \cJ_m\,,\quad f_i\in R\,,
\end{equation}

where $f_i$ are the ring elements, i.e., they are polynomials of the parameters. However, since $\cJ_s$ and $\cJ_{1,\dots,m}$ are elements in the invariant ring, $f_i$ should also be invariants. They can be parameterized by the generating invariants in a polynomial form

\begin{equation}
\label{eq:jpower_new}
\cJ_s^t=\sum_{i=1}^{m}P_i(\cJ_l,\cI_k)\cJ_i\,,\quad i\neq s\,.
\end{equation}

Therefore, Hilbert's Nullstellensatz tells us that if a CP-odd invariant $\cJ_s$ is redundant in the presence of a given CPC set, there must exist a syzygy of some power of $\cJ_s$ that can be used to eliminate this invariant. This theorem is quite helpful when using the syzygy approach. For example, in Ref. <cit.>, the syzygies are used to identify the CPC conditions in the 2HDM, one can find that the syzygies regarding the vanishing invariant $\cJ_s$ can be obtained at order $\cJ_s^2$, there is no need to discuss the spectrum of the syzygies as presented in their analysis.

It is possible to come up with an elimination algorithm based on the Hilbert's Nullstellensatz to find the minimal CPC set. However, the problem of finding a syzygy like Eq. (<ref>), or in more mathematical language, determining whether an ideal $\cJ_s^t$ is a subset of another ideal $I:=(\cJ_1,\dots,\cJ_m)$ highly relies on the calculation of the Gröbner basis, which is computationally quite expensive in complicated polynomial rings. The undetermined power $t$ also introduces a lot of complexity in this problem. There are different methods and packages that are suitable for this problem. As outlined in App. <ref>, the syzygy problem can also be converted to a finite system of linear equations, which can be addressed using standard linear algebra techniques. However, some high-degree syzygies lead to very complex linear systems that can not be easily solved. Some software systems are devoted to studying the algebraic geometry and commutative algebra, such as  <cit.> and  <cit.>. They have proved helpful when exploring the algebraic structures of simple theories. In addition, there is the  function PolynomialReduce, that can also be used to solve these problems. However, it is also based on the expensive Gröbner basis calculation. Despite the capabilities of these packages, they failed to generate results within a reasonable time frame for the three-generation case of the .

Although we were unable to construct the minimal CPC set for the theory with three generations of fermions, we did identify some example of syzygy that follows the form shown in Eq. (<ref>). For instance,

\begin{equation}
\begin{split}
2 \cJ_{14}^2 &= \cJ_1 \left(2 \cJ_{34} \cI_2-2 \cJ_{76}-\cJ_{11} \cI_2^2-\cJ_{11} \cI_5\right) + \cJ_3 \left(4 \cJ_{36}- 4 \cJ_{13} \cI_2 + \cI_2^2\cJ_3+\cI_5\cJ_3\right)\\
&+2\cJ_4 \left(\cJ_{11} \cI_2 - \cJ_{34}\right)-2 \cJ_{11} \cJ_{15} + 2 \cJ_{13}^2\,,
\end{split}
\end{equation}

one can find that if $\{\cJ_1,\cJ_3,\cJ_4,\cJ_{11},\cJ_{13}\}$ is set to zero, then $\cJ_{14}=0$ automatically, which means if we include the list of five invariants in the CPC set, the inclusion of $\cJ_{14}$ becomes unnecessary.

In addition, we also find that some CP-odd invariants must be added to the minimal CPC set. Specifically, if only one CP-odd invariant $\cJ_s$ in the generating set is non-zero under a specific spectrum, it would not be possible to establish Eq. (<ref>). Consequently, $\cJ_s$ cannot be eliminated by any CPC set unless it is included in the set. We have found 3 such cases where only one CP-odd invariant is non-zero. The first is the limit $M_N \to 0$ where the theory is reduced to a copy of the SM quark sector. Then the analogue of the Jarlskog invariant $\mathcal{J}_{10}=\,\Im\Tr(\neat{X_{N}^{2} X_{e}^{2} X_{N}^{} X_{e}^{}})$ is the only non-vanishing CP-odd invariant and has to be included in the minimal CPC set. The other 2 cases are found under the spectrum $\{M_N \to m_N \mathbb{1}, Y_e \to 0\}$ and $\{M_N \to m_N \mathbb{1},Y_N \to y_N \mathbb{1}\}$ which force us to add $\mathcal{J}_{74}=\,\Im\Tr(\neat{Y_{N}^{} Y_{N}^{\dagger} Y_{N}^{} M_{N}^{*} Y_{N}^{T} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger} Y_{N}^{} M_{N}^{*} Y_{N}^{T} Y_{N}^{*} Y_{N}^{T} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}})$ and $\mathcal{J}_{251}=\,\Im\Tr(\neat{X_{e}^{2} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*2} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger} X_{e}^{} Y_{N}^{} M_{N}^{*} Y_{N}^{T} X_{e}^{*} Y_{N}^{*} M_{N}^{} Y_{N}^{\dagger}})$ to the minimal CPC set respectively. Here, as before, we define $X_{e,N} = \neat{Y_{e,N}^{}Y_{e,N}^{\dagger}}$. It's interesting that we have to rely on the highest-order invariant $\mathcal{J}_{251}$ to establish the CPC conditions of the theory. Notably, $\mathcal{J}_{251}$ is the last invariant in the CP-odd generating set, as shown in App. <ref>.

We want to emphasize one point here. The parameters of the neutrino sector are not measured well enough today to even exclude that one of the neutrinos is massless. Therefore, having a set of flavor invariants which determines the CPC conditions without assumptions on the spectrum is important to make general statements about the theory that hold true for all possible experimental results.

In conclusion, we introduced Hilbert's Nullstellensatz as a reliable mathematical language to address the problem of the CPC condition. However, the problem is inherently complex and cannot be solved with our current efforts. It deserves further exploration as it requires larger computational resources to be successfully resolved. Additionally, we have solved the complete CPC conditions of the for the two-generation case, as detailed in App. <ref>, and some special conditions corresponding to the pseudo-real couplings are found, which deserve further investigation.

§ CONCLUSIONS

In this paper, we have analyzed the algebraic structure of the , the SM extended with three generations of sterile neutrinos, by constructing its set of generating flavor invariants. In the quark sector, the theory exhibits a complete intersection ring, and the flavor invariants can be easily generated with the guidance of the Hilbert series. However, in the , the theory corresponds to a non-complete intersection ring, introducing many complexities. The analysis of the Hilbert series in this context requires a more careful examination to extract correct information about the number of generating invariants and syzygies. We introduce a novel graph-based approach for constructing flavor invariants. This method maps all single-trace invariants to walks on the graph, allowing for a systematic construction of invariants. Notably, this method can be extended to construct flavor invariants (or covariants) in broader theoretical frameworks beyond the scope of our current work. The resulting invariants are then reduced to a generating set by removing redundant invariants. The number of generating invariants and syzygies are carefully compared with the plethystic logarithm and we find agreement at each order in the flavor spurions. The generating set comprises of in total 459 invariants, out of which 208 are CP-even and 251 are CP-odd. We have furthermore reduced this set to a primary set of flavor invariants capturing the 21 physical parameters of the renormalizable Lagrangian.

We can show that the generating set is able to capture the physical degrees of freedom of the seesaw limit. Considering structures in the flavor invariants that are suppressed with one and two powers of the heavy Majorana mass in the seesaw limit, we showed that the 12 and 21 physical parameters are captured by our generating set of invariants at the respective orders. Here, the number 12 accounts for the number of masses, mixing angles and phases in the low-energy theory of neutrinos and charged leptons. By also taking the total suppression of the invariants with the heavy Majorana mass in the effective theory into account, it can be seen in a straightforward way, that the 12 and 21 physical parameters can only be accessed with a suppression of two and four powers of the Majorana mass respectively.

In the spirit of the Jarlskog invariant, we have defined a minimal set of CP-odd invariants that can determine the CPC conditions of the theory, which is called the minimal CPC set. This set is defined to be assumption-free regarding the spectrum of the parameters, making it available to determine whether CPV exists in the theory for all upcoming experimental results. The reduction of the CP-odd generating set to a minimal CPC set turns out to be a hard task and we only managed to reduce the set of CP-odd flavor invariants to a more minimal set in the case of two generations of fermions, while the reduction of the three generation case was too complicated to complete by our means. Along the way, we present useful algorithms which are essential to our analysis and hopefully prove to be useful for other analyses of flavor invariants in the future. It is worth noting that the flavor structures of the type-I and type-III seesaws are identical to that of the , while the couplings of the type-II seesaw cannot be identified with those of the due to the different particle content. Instead, the only flavorful parameter in the theory has the same transformation properties as those studied in Ref. <cit.>. Therefore, the flavor invariants identified in this paper can be directly applied to the type-I and type-III seesaw models by renaming the flavor matrices.

An obvious next step for the flavor invariants is to connect them to phenomenological applications to study, for instance, CPV observables in a consistent way with the help of flavor invariants.

§ ACKNOWLEDGMENTS

We thank Quentin Bonnefoy, Gabriele Dian, Emanuele Gendy, Gabriel Massoni Salla, Paula Pilatus and Pablo Quílez for useful discussions. DL thanks the DESY theory group for their hospitality during a long-term visit, where part of this work was completed. He also thanks the ENS de Lyon for providing the scholarship that made this visit possible. This work is supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy EXC 2121 “Quantum Universe” – 390833306, as well as by the grant 491245950. This project also has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie Staff Exchange grant agreement No 101086085 - ASYMMETRY. C.Y.Y. is supported in part by the Grants No. NSFC-11975130, No. NSFC-12035008, No. NSFC-12047533, the National Key Research and Development Program of China under Grant No. 2017YFA0402200, the China Postdoctoral Science Foundation under Grant No. 2018M641621 and the Helmholtz-OCPC International Postdoctoral Exchange Fellowship Program.

§ APPENDICES

In App. <ref>, we present useful parameterizations for the flavorful couplings that we use throughout this paper. In App. <ref>, we show our results for the multi-graded Hilbert series and the multi-graded plethystic logarithm that we do not show in the main text due to their length. In App. <ref>, we present the algorithms used in this paper, which include the algorithm to generate flavor invariants and the algorithm to reduce the invariants to a generating set. In App. <ref>, we list all invariants in the generating set along with their CP parities. In App. <ref>, we introduce the Hironaka decomposition that we use in some parts of the paper. In App. <ref>, we introduce Hilbert's Nullstellensatz (a fundamental theorem in algebraic geometry stating the relationship between polynomial equations and the points where they vanish) and present the relevant concepts that are useful for obtaining a minimal set of CP-odd invariants that capture the necessary and sufficient conditions for CPC. Finally, in App. <ref>, we present the complete list of CPC conditions for the two generation case obtained through the ideal-related method.

§ PARAMETERIZATION OF FLAVOR MATRICES

§.§ Standard parameterization

The parameterizations of the flavor matrices in the lepton and quark sectors have been discussed in detail in Refs. <cit.>. We use a slightly different convention in this paper and thus, we will follow a similar discussion here to introduce the parameterization in this new convention. The fermionic part of the Lagrangian is given by

\begin{equation}
\mathcal{L}_4
\supset\sum_{\Psi}\bar{\Psi}i \slashed{D}\Psi
-\left[\frac{1}{2}\left(NCM_NN\right)+\bar{L}Y_NN\widetilde{H}+\bar{L}Y_e e H
+\bar{Q}Y_u u \widetilde{H}+\bar{Q}Y_d d H +\mbox{H.c.}\right],
\end{equation}

where the kinetic term sums all fermion fields $\Psi=\{Q, L, u, d, e, N\}$.
The flavor transformations of the fermion fields are given by

\begin{equation}
\begin{split}
L& \to U_L\, L\,, \quad  e \to U_e\, e\,, \quad  N \to U_N\, N\,, \\
Q& \to U_Q\, Q\,, \quad  u \to U_u\, u\,, \quad  d \to U_d\, d\,. \\
\end{split}
\end{equation}

The corresponding Yukawa matrices and the Majorana mass matrix transform as

\begin{equation}
\label{eq:flavor_transformation}
\begin{split}
Y_e& \to U_L Y_e U_e^{\dagger}\,, \quad Y_N \to U_L Y_N U_N^{\dagger}\,,\quad M_N \to U_N^* M_N U_N^{\dagger}\,, \\
Y_u& \to U_Q Y_u U_u^{\dagger}\,, \quad Y_d \to U_Q Y_d U_d^{\dagger} \,.\\
\end{split}
\end{equation}

The Yukawa matrices $Y_{N,e,u,d}$ are general complex matrices, while the Majorana mass matrix $M_N$ is symmetric, which can be diagonalized as follows

\begin{equation}
\label{eq:diagonalization}
\begin{split}
\end{split}
\end{equation}

where $\widehat{Y}_{N,e,u,d}$ and $\widehat{M}_N$ are diagonal matrices with real and non-negative entries. $V_{f},W_{f}$ with $f=N,e,u,d$ and $V'_N$ are unitary matrices. We can choose specific flavor transformations in Eq. (<ref>) to get a fixed mass basis. In this paper, we will work on the charged lepton diagonal basis in the lepton sector and up basis in the quark sector. This can be achieved by setting $U_L=V_e^{\dagger}, U_e=W_e^{\dagger}, U_N=V_N'^{T}, U_Q=V_u^{\dagger}, U_u=W_u^{\dagger}, U_d=W_d^{\dagger}$. The flavor matrices will be fixed to the following forms,

\begin{equation}
\label{eq:basis}
\begin{split}
\end{split}
\end{equation}

where $V_\CKM=V_u^{\dagger}V_d$ is the CKM matrix, corresponding to the mismatch between the diagonalization matrices of the up and down sectors. $V_{L}=V_e^{\dagger}V_N$ is a similar matrix in the lepton sector, describing the mismatch between the diagonalization matrices of lepton Yukawa matrices. The existence of the Majorana mass matrix in the lepton sector introduces another mixing matrix $W=V_N'^{T}W_N$, which describes the mismatch between the diagonalization matrices of $Y_N$ and $M_N$. If the Majorana mass term is forbidden in the theory, the neutrino will have a Dirac mass and we can simply choose $U_N=W_N$, which will lead to $W=\mathbbm{1}$. The lepton sector will have the same structure as the quark sector and $V_L$ can be identified with the PMNS matrix, phenomenologically.

From Eq. (<ref>), if we assume that the flavor matrices have no degenerate or vanishing eigenvalues,[If there is a degenerate or vanishing mass spectrum, there will be a larger degrees of freedom to redefine the mixing matrix, some parameters in our current parameterization could be unphysical, and they can be removed by those redefinition. Under this enlarged symmetry, we should adopt a new parameterization, the number of parameters should exactly match the number of physical observables. On the other hand, if we adopt the most general parameterization, and set the special spectrum afterwards, some parameters in the parameterization will be redundant. In some cases, the redundant parameters could be canceled in the expression, and they never appear in the flavor invariants. In other cases, some parameters will appear together as a single polynomial in different invariants. In the latter case, the polynomial that correlates them should be considered as a single parameter in the parameterization, which should remove the redundancies.] the diagonalization matrices will at least have a column phase redefinition freedom,  $V_{f}\rightarrow V_f e^{i\widehat{\Phi}_f}\,, W_{f}\rightarrow W_f e^{i\widehat{\Phi}_f}, f=\{N,e,u,d\}$ and $V_N'\rightarrow V_N'\eta_N$. This can be expressed as the a rephasing invariance of the diagonal matrices

\begin{equation}
\begin{aligned}
e^{i \widehat{\Phi}_f} \widehat{Y}_f e^{-i \widehat{\Phi}_f} & = \widehat{Y}_f, \quad f=N,e,u,d\,, \\
\eta_N \widehat{M}_N \eta_N & = \widehat{M}_N \,,
\end{aligned}
\end{equation}

where $\widehat{\Phi}_{f}\equiv \diag(\phi_{f_1}, \phi_{f_2}, \phi_{f_3})\,, f=N,e,u,d$ are diagonal complex phase matrices, and $\eta_N$ is a diagonal matrix with eigenvalues $\pm 1$. Under these rephasings, the mixing matrices transform as

\begin{equation}
\label{eq:phase_redefinition}
V_\CKM \rightarrow e^{-i \widehat{\Phi}_u}V_\CKM e^{i\widehat{\Phi}_d}\,,\quad
V_L \rightarrow e^{-i \widehat{\Phi}_e}V_L e^{i \widehat{\Phi}_N}\,,\quad
W \rightarrow  \eta_N W e^{i \widehat{\Phi}_N}\,,
\end{equation}

The $3\times 3$ unitary matrix can be parameterized as

\begin{equation}
\label{eq:unitary_matrix}
U_3= e^{i \varphi}e^{i \widehat{\Psi}}U(\theta_{12},\theta_{13},\theta_{22},\delta)e^{i \widehat{\Phi}}\,,
\end{equation}

where $\varphi$ is a overall phase, $\widehat{\Psi}=\diag(0,\psi_1,\psi_2)$ and $\widehat{\Phi}=\diag(0,\phi_1,\phi_2)$, and $U(\theta_{12},\theta_{13},\theta_{23},\delta)$ takes the standard form

\begin{eqnarray}
\label{eq:ckm}
&& U(\theta_{12},\theta_{13},\theta_{23},\delta) \equiv
\left( \begin{array}{ccc}
1 & 0 & 0 \\
0 & c_{23} & s_{23} \\
0 & -s_{23} & c_{23} \end{array} \right)
\left( \begin{array}{ccc}
c_{13} & 0 & s_{13}e^{-i \delta} \\
0 & 1 & 0 \\
- s_{13}e^{i \delta} & 0 & c_{13} \end{array} \right)
\left( \begin{array}{ccc}
c_{12} & s_{12} & 0 \\
-s_{12} & c_{12} & 0\\
0 & 0 & 1 \end{array} \right)\,,
\end{eqnarray}

where $s_{ij}\equiv\sin\theta_{ij}$, $c_{ij}\equiv\cos\theta_{ij}$, and $\theta_{ij}\in [0,\pi/2], \delta\in [0,2\pi)$.

In the quark sector, by introducing the phase redefinition in Eq. (<ref>), the phases in the unitary matrix $V_\CKM$ can be  absorbed, it takes the standard form as given in Eq. (<ref>). The number of parameters are summarized as follows

\begin{eqnarray*}
\begin{array}{c|ccc}
\hline
\text{Matrices}& \text{Masses} & \text{Angles} & \text{Phases} \\
\hline
\widehat{Y}_u &  3 & 0 & 0  \\
\widehat{Y}_d &  3 & 0 & 0  \\
V_{\CKM} & 0 & 3 & 1\\[5pt]\hline
\text{Total} & 6 &  3 & 1\\
\hline
\end{array}
\end{eqnarray*}

There are a total of 10 parameters, consisting of 6 quark masses, 3 mixing angles, and 1 CP phase.

In the lepton sector, we can parameterize $V_L$ and $W$ as

\begin{equation}
\begin{split}
\label{eq:VW}
V_L=&\, e^{i \varphi}e^{i \widehat{\Psi}}U(\theta_{12},\theta_{13},\theta_{23},\delta)e^{i \widehat{\Phi}}\,,\\
W=&\, e^{i \varphi'}e^{i \widehat{\Psi}'/2}U(\theta_{12}',\theta_{13}',\theta_{23}',\delta')e^{i \widehat{\Phi}'}\,.
\end{split}
\end{equation}

Since $V_L$ and $W$ share the same rephasing matrix $\widehat{\Phi}_N$, we can use this freedom to remove either $\widehat{\Phi}$ or $\widehat{\Phi}'$ in Eq. (<ref>), depending on what observable we are interested in. If the flavor invariant only depends on $Y_N$ and $Y_e$, then the mixing matrix $W$ can be simply set to $W=\mathbbm{1}$, since it is not an observable. We can use the phase matrices $\widehat{\Phi}_e$ and $\widehat{\Phi}_N$ to remove the phases of the unitary matrix $V_L$, then it takes the same form as given in Eq. (<ref>). The mixing matrices in this case take the following form

\begin{equation}
V_L=U(\theta_{12},\theta_{13},\theta_{23},\delta), \quad W=\mathbbm{1}\,.
\end{equation}

If the flavor invariant only depends on $Y_N$ and $M_N$, then the mixing matrix $V_L$ is not an observable, which can be set to identity matrix. We can use the rephasing matrix $\widehat{\Phi}_N$ to remove the phases $\varphi'$ and $\widehat{\Phi}'$ in Eq. (<ref>), the effect of $\eta_N$ is indicated by the factor of two in the phase matrix $e^{i\widehat{\Psi}'/2}$, which limits the phases in the range of $\psi'_i/2\in [0,\pi)$ if by convention $\psi'_i\in [0,2\pi)$. In this case the mixing matrices will take the following form

\begin{equation}
V_L = \mathbbm{1}\,,\quad W=\, e^{i \widehat{\Psi}'/2}U(\theta_{12}',\theta_{13}',\theta_{23}',\delta')\,.
\end{equation}

If the flavor invariant depends on the three flavor matrices $Y_e, Y_N$ and $M_N$,[The flavor invariants can not be constructed by only $Y_e$ and $M_N$, since they are disconnected objects if $Y_N$ is not included, as shown in Fig. <ref>.] we can use the rephasing phases $\widehat{\Phi}_e$ and $\widehat{\Phi}_N$ to remove the phases in $V_L$ and the phase $\varphi'$ in $W$ in Eq. (<ref>), without loss of generality. Then, the mixing matrices take the following form

\begin{equation}
V_L = U(\theta_{12},\theta_{13},\theta_{23},\delta)\,, \quad W=\, e^{i \widehat{\Psi}'/2}U(\theta_{12}',\theta_{13}',\theta_{23}',\delta')e^{i \widehat{\Phi}'}\,.\quad
\end{equation}

We can also use $\widehat{\Phi}_e$ and $\widehat{\Phi}_N$ to remove the phases $e^{i\varphi}e^{i\widehat{\Psi}}$ and $e^{i\varphi'}e^{i\widehat{\Phi}'}$ in $V_L$ and $W$ in Eq. (<ref>), respectively, then the mixing matrices are given by

\begin{equation}
\label{eq:param_lepton}
V_L = U(\theta_{12},\theta_{13},\theta_{23},\delta)e^{i\widehat{\Phi}}\,, \quad W=\, e^{i \widehat{\Psi}'/2}U(\theta_{12}',\theta_{13}',\theta_{23}',\delta')\,.\quad
\end{equation}

The third case corresponds to the most general parameterization of the mixing matrices $V_L$ and $W$ and we can work in this basis without loss of generality. To be specific, we use the following parameterization

\begin{equation}
\label{eq:parameterization_3gen}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\ Y_N = V_L \cdot \diag\(y_{1},y_{2},y_{3}\) \cdot W^{\dagger},\  M_N = \diag\(m_1,m_2,m_3\) \, ,
\end{equation}

with $V_L$ and $W$ defined in Eq. (<ref>). We summarize the the number of parameters in each matrix in Tab. <ref>.

For a more comprehensive discussion, please refer to Ref. <cit.>. In their paper, the cases for $n_N=n_f=2$ and $n_N=2,n_f=3$ in the lepton sector are also discussed. For $n_N=n_f=2$, all flavor matrices are $2\times 2$ matrices, the $2\times2$ unitary matrix is given by

\begin{equation}
U_2=e^{i\varphi} \diag(1,e^{i\psi})\cdot
\begin{pmatrix}
\cos\theta & \sin\theta\\
-\sin\theta & \cos\theta
\end{pmatrix}\cdot\diag(1,e^{i\phi})\,.
\end{equation}

Follow a similar discussion, we can easily find the following parameterization,

\begin{equation}
\label{eq:parameterization_2gen}
Y_e = \diag\(y_e,y_{\mu}\), \qquad Y_N = V_L \cdot \diag\(y_{1},y_{2}\) \cdot W^{\dagger}, \qquad M_N = \diag\(m_1,m_2\) \, ,
\end{equation}

\begin{equation}
\label{eq:mixing_matrix_2gen}
V_L =\begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \cdot \diag\(1,e^{i \phi}\), \quad W = \diag\(1,e^{i \varphi}\) \cdot \begin{pmatrix} \cos\alpha & \sin\alpha \\ -\sin\alpha & \cos\alpha \end{pmatrix} \, .
\end{equation}

For $n_N=2,n_f=3$, $Y_N$ is a $3\times 2$ matrix, and $M_N$ is a $2\times 2$ matrix. Similarly, the flavor matrices can be parameterized as

\begin{equation}
\label{eq:parameterization_3x2gen}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N = V_L \cdot
\begin{pmatrix}
y_1 & 0\\
0 & y_2\\
0 & 0
\end{pmatrix} \cdot W^{\dagger},\quad  M_N = \diag\(m_1,m_2\) \, ,
\end{equation}

with $V_L$ and $W$ defined as follows

\begin{equation}
\label{eq:mixing_matrix_3x2gen}
V_L =U(\theta_{12},\theta_{13},\theta_{23},\delta)\cdot \diag \(1,e^{i \phi},1\)\,, \quad
W = \diag\(1,e^{i \varphi}\) \cdot \begin{pmatrix} \cos\alpha & \sin\alpha \\ -\sin\alpha & \cos\alpha \end{pmatrix} \,,
\end{equation}

which is analogous to the parameterization in Eq. (<ref>), however the last phase in $\hat{\Phi}$ is unphysical, since it vanishes when multiplied with the zeros in the last row of the $3\times2$ block diagonal matrix.

\begin{eqnarray*}
\begin{array}{c|ccc}
\hline
\text{Matrices}& \text{Masses} & \text{Angles} & \text{Phases} \\ \hline
\widehat{M}_N &  3\,[2]\,(2) & 0 & 0  \\
\widehat{Y}_N &  3\,[2]\,(2) & 0 & 0  \\
\widehat{Y}_e &  3\,[2]\,(3) & 0 & 0  \\
V_L & 0 & 3\,[1]\,(3) & 3\,[1]\,(2) \\
W & 0 & 3\,[1]\,(1) & 3\,[1]\,(1)  \\
\hline
\text{Total} & 9\,[6]\,(7) &  6\,[2]\,(4) & 6\,[2]\,(3)\\\hline
\end{array}
\end{eqnarray*}

The number of masses, mixing angles and phases in the parameterization of the lepton sector for the case of $n_N=n_f=3\,[n_N=n_f=2]\,(n_N=2,n_f=3)$.

In the Dirac limit, $M_N$ is set to zero, and $W$ is not relevant, we can simply set $W=\mathbbm{1}$. In the case of $n_N=n_f=3$, the matrix $V_L$ will become the PMNS matrix $V_{\PMNS}$, which takes the same standard form as $V_{\CKM}$. The parameterization is given as

\begin{equation}
\label{eq:parameterization_3gen_Dirac1}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N = U(\theta_{12},\theta_{13},\theta_{23},\delta) \cdot \diag\(y_{1},y_{2},y_{3}\)\, .
\end{equation}

Similarly, for $n_N=n_f=2$, it can be parameterized as

\begin{equation}
\label{eq:parameterization_2gen_Dirac1}
Y_e = \diag\(y_e,y_{\mu}\), \qquad Y_N = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \cdot \diag\(y_{1},y_{2}\)\, ,
\end{equation}

where the phase of $V_L$ in the Majorana case can be removed by the rephasing of the right-handed neutrino $N$. For $n_N=2, n_f=3$, the parameterization is given by

\begin{equation}
\label{eq:parameterization_3x2gen_Dirac1}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N = U(\theta_{12},\theta_{13},\theta_{23},\delta) \cdot
\begin{pmatrix}
y_1 & 0\\
0 & y_2\\
0 & 0
\end{pmatrix}\, ,
\end{equation}

where the phase $e^{i\phi}$ in Eq. (<ref>) can be absorbed by the field $N$ in the Dirac case. The parameters in the Dirac limit are summarized in Tab. <ref>.

\begin{eqnarray*}
\begin{array}{c|ccc}
\hline
\text{Matrices}& \text{Masses} & \text{Angles} & \text{Phases} \\
\hline
\widehat{Y}_N &  3[2](2) & 0 & 0  \\
\widehat{Y}_e &  3[2](3) & 0 & 0  \\
V_{\PMNS}  & 0 & 3[1](3) &1[0](1) \\
\hline
\text{Total} & 6[4](5) & 3[1](3) & 1[0](1)\\\hline
\end{array}
\end{eqnarray*}

The number of masses, mixing angles and phases in the parameterization of the lepton sector for the case of $n_N=n_f=3\,[n_N=n_f=2]\,(n_N=2,n_f=3)$ in the Dirac limit.

§.§ Algebraic parameterization

The parameterization described above is favored for its phenomenological relevance, as Yukawa matrices are factorized into eigenvalues and mixing matrices, which aligns with experimental observables. However, when exploring the algebraic structures of invariants, the inclusion of trigonometric functions introduces complexity. The use of sine and cosine as objects in the polynomial expansion of flavor invariants can complicate the exploration of these structures.  Consequently, alternative parameterizations that are more suitable for polynomial expressions are needed.

One possible solution is to parameterize the trigonometric functions. A frequently used parameterization for the unit circle is provided as follows:

\begin{equation}
x(t) = \frac{1-t^2}{1+t^2},\quad y(t) = \frac{2t}{1+t^2}\,,\quad \text{with~} t\in(-\infty,+\infty).
\end{equation}

However the point $(-1,0)$ on the unit circle can only be obtained in the limit $t\to\infty$. Another parameterization that can cover the whole circle is given by

\begin{equation}
x(t) = \frac{1-6t^2+t^4}{1+2t^2+t^4},\quad y(t) = \frac{4t-4t^3}{1+2t^2+t^4}\,,\quad \text{with~} t\in (-1,1]\,.
\end{equation}

Thus the sine and cosine functions in the mixing matrices of our above parameterization can be replaced with $x(t)$ and $y(t)$ respectively. However, this parameterization may introduce new complexities, as it leads to rational polynomials.

Following the parameterization in Eq. (<ref>), it is possible to work with the diagonal basis of both $Y_e$ and $M_N$ while leaving $Y_N$ undiagonalized. The parameterization is given as follows

\begin{equation}
\label{eq:parameterization_3gen_new}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N = \begin{pmatrix}
r_{11} & c_{12} & c_{13}\\
r_{21} & c_{22} & c_{23}\\
r_{31} & c_{32} & c_{33}\\
\end{pmatrix},\quad
M_N = \diag\(m_1,m_2,m_3\) \, ,
\end{equation}

where $r_{ij}$ and $c_{ij}$ are general labels for real and complex parameters respectively, and $c_{ij}$ can be easily written as real and imaginary parts by introducing two real parameters. The phases of the first column of $Y_N$ are absorbed by the rephasing of charged lepton fields.[It is not necessary for the rephasing degree of freedom to target the first column, any phase in each row of $Y_N$ can be eliminated.] In this parameterization all invariants are expressed as polynomials of simple variables, making it easier to analyze the algebraic structures of the theory. It is easy to see that the number of parameters in this parameterization is still 21, which is the same as the physical parameterization in Eq. (<ref>). Similarly, for $n_N=n_f=2$, we can parameterize it as

\begin{equation}
\label{eq:parameterization_2gen_new}
Y_e = \diag\(y_e,y_{\mu}\), \qquad Y_N = \begin{pmatrix}
r_{11} & c_{12}\\
r_{21} & c_{22}
\end{pmatrix}, \qquad
M_N = \diag\(m_1,m_2\) \, .
\end{equation}

For $n_N=2, n_f=3$, the parameterization is given by

\begin{equation}
\label{eq:parameterization_3x2gen_new}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N =
\begin{pmatrix}
r_{11} & c_{12}\\
r_{21} & c_{22}\\
r_{31} & c_{32}
\end{pmatrix},\quad  M_N = \diag\(m_1,m_2\) \, .
\end{equation}

The mapping between these two parameterizations can be found by solving equations built from the entries of $Y_N$, which will not be shown here.

Under the Dirac limit, the Majorana mass term $M_N$ is set to 0, the parameters in Yukawa matrix $Y_N$ can be further reduced by the field redefinition of right-handed neutrino $N$. Starting from the $Y_N$ in Eq. (<ref>), we can absorb the phases in $c_{12}$ and $c_{13}$ by rephasing of $N$, these two real parameters can be further set to 0 by two rotations of $N$, i.e., $N\to R_{12}(\theta_2)R_{23}(\theta_1)N$, where $R_{ij}$ is the rotation matrix acting on the $(i,j)$ entries of $N$. The rotation $R_{12}$ can make $r_{21}$ and $r_{31}$ complex, but their phases can always be absorbed by the rephasing of charged lepton fields. With these two leading zeros in the second and third columns, the rephasing of $N$ can be used again to remove the phases in $c_{22}$ and $c_{23}$, the corresponding real parameters $r_{22}$ and $r_{23}$ can be mixed by the rotation $N\to R_{23}(\theta_3)N$. By setting proper value of $\theta_3$, the parameter $r_{23}$ can be set to zero. The leading two zeros in the third column making it possible to remove the phase in $c_{33}$ by rephasing of $N$, resulting in a real parameter $r_{33}$. The whole transformation process can be presented as follows

\begin{equation}
Y_N\to \begin{pmatrix}
r_{11} & r_{12} & r_{13} \\
r_{21} & c_{22} & c_{23} \\
r_{31} & c_{32} & c_{33} \\
\end{pmatrix}\to
\begin{pmatrix}
r_{11} & 0 & 0 \\
r_{21} & c_{22} & c_{23} \\
r_{31} & c_{32} & c_{33} \\
\end{pmatrix} \to
\begin{pmatrix}
r_{11} & 0 & 0 \\
r_{21} & r_{22} & r_{23} \\
r_{31} & c_{32} & c_{33} \\
\end{pmatrix}\to
\begin{pmatrix}
r_{11} & 0 & 0 \\
r_{21} & r_{22} & 0 \\
r_{31} & c_{32} & c_{33} \\
\end{pmatrix}\to
\begin{pmatrix}
r_{11} & 0 & 0 \\
r_{21} & r_{22} & 0 \\
r_{31} & c_{32} & r_{33} \\
\end{pmatrix}\,.
\end{equation}

There is no further field redefinition that can be used to reduce the number of parameters, as a result, in the Dirac limit with $n_N=n_f=3$, we can parameterize the flavor matrices as

\begin{equation}
\label{eq:parameterization_3gen_Dirac}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N = \begin{pmatrix}
r_{11} & 0 & 0\\
r_{21} & r_{22} & 0\\
r_{31} & c_{32} & r_{33}\\
\end{pmatrix}\, .
\end{equation}

We find there are exactly 10 real parameters in this parameterization, as expected. Similarly, for $n_N=n_f=2$ the flavor matrices are parameterized as

\begin{equation}
\label{eq:parameterization_2gen_Dirac}
Y_e = \diag\(y_e,y_{\mu}\), \qquad Y_N = \begin{pmatrix}
r_{11} & 0\\
r_{21} & r_{22}
\end{pmatrix}\, .
\end{equation}

For $n_N=2, n_f=3$ the parameterization is given by

\begin{equation}
\label{eq:parameterization_3x2gen_Dirac}
Y_e = \diag\(y_e,y_{\mu},y_{\tau}\),\quad Y_N =
\begin{pmatrix}
r_{11} & 0\\
r_{21} & r_{22}\\
r_{31} & c_{32}
\end{pmatrix}\, .
\end{equation}

§ RESULTS FOR MULTI-GRADED HILBERT SERIES AND PLETHYSTIC LOGARITHM

§.§ Model with nN=nf=3

In Eqs. (<ref>) and (<ref>), we have presented the ungraded Hilbert series, where the same grading is used for all spurions. However, the information encoded in the ungraded HS is not enough for analyses, where the PL is used to count the number of generating invariants and syzygies at each order. Therefore, in this section, we present the multi-graded HS, where the single spurion $q$ in Eqs. (<ref>) and (<ref>) is traded for multiple spurions $e,m$ and $n$, corresponding to the spurions of the flavor matrices $Y_e,M_N$ and $Y_N$ respectively. The denominator of the multi-graded HS is given as follows

\begin{equation}\label{eq:gradedHS_den}
\begin{split}
\cD&(e,m,n)=(1-e^2)(1-m^2)(1-n^2)(1-e^4)(1-m^4)(1-n^4)(1-e^2 n^2)^2(1-m^2 n^2)\\
&(1-e^6)(1-m^6)(1-n^6)(1-e^2 n^4)(1-e^4 n^2)(1-m^2 n^4)^2(1-m^4 n^2)(1-e^2 m^2 n^2)\\
&(1-m^2 n^6)(1-m^4 n^4)(1-e^2 m^2 n^4)(1-e^2 m^4 n^2)(1-e^4 m^2 n^2)(1-m^2 n^8)\\
&(1-e^2 m^2 n^6)(1-e^4 m^2 n^4)(1-e^4 m^4 n^2)(1-e^4 m^2 n^6)(1-e^4 m^4 n^4)(1-e^6 m^2 n^4)\\
&(1-e^4 m^2 n^8)(1-e^8 m^2 n^4)(1-e^8 m^4 n^4)\,,
\end{split}
\end{equation}

while there are 6582 terms in the numerator, which goes up to order $\cO([emn]^{196})$ in total powers of spurions. Due to its length we only show the terms up to $\cO([emn]^{26})$ below

\begin{align*}
\cN&(e,m,n) = 1 - e^2 n^2 + 2 m^4 n^4 + e^4 n^4 + 2 e^2 m^2 n^4 + 2 m^4 n^6 + 2 m^6 n^4 + 4 e^2 m^2 n^6\\
&+ 4 e^2 m^4 n^4 + 3 e^4 m^2 n^4 + 3 m^4 n^8 + 3 m^6 n^6 + m^8 n^4 + 5 e^2 m^2 n^8 + 7 e^2 m^4 n^6 + 3 e^2 m^6 n^4\\
&+ 4 e^4 m^2 n^6 + 5 e^4 m^4 n^4 + e^6 m^2 n^4 + m^4 n^{10} + 3 m^6 n^8 + m^8 n^6 + 3 e^2 m^2 n^{10} + 9 e^2 m^4 n^8\\
& + 6 e^2 m^6 n^6 + e^2 m^8 n^4 + 4 e^4 m^2 n^8 + 10 e^4 m^4 n^6 + 5 e^4 m^6 n^4 + 3 e^6 m^2 n^6 + 3 e^6 m^4 n^4\\
& + m^4 n^{12} + m^6 n^{10} + 3 m^8 n^8 + e^2 m^2 n^{12} + 5 e^2 m^4 n^{10} + 9 e^2 m^6 n^8 + 2 e^2 m^8 n^6 + e^4 m^2 n^{10}\\
& + 16 e^4 m^4 n^8 + 10 e^4 m^6 n^6 + 2 e^4 m^8 n^4 + 4 e^6 m^2 n^8 + 8 e^6 m^4 n^6 + 3 e^6 m^6 n^4 + 2 e^8 m^2 n^6\\
& + e^8 m^4 n^4 + m^8 n^{10} + m^{10} n^8 + e^2 m^4 n^{12} + 6 e^2 m^6 n^{10} + 7 e^2 m^8 n^8 - e^2 m^{10} n^6\numberthis\\
& + 10 e^4 m^4 n^{10} + 19 e^4 m^6 n^8 + 4 e^4 m^8 n^6 + e^6 m^2 n^{10} + 14 e^6 m^4 n^8 + 8 e^6 m^6 n^6 + e^6 m^8 n^4\\
& + e^8 m^2 n^8 + 5 e^8 m^4 n^6 + 2 e^8 m^6 n^4 + e^{10} m^2 n^6 - m^6 n^{14} + 2 m^{10} n^{10} + m^{12} n^8 - 3 e^2 m^4 n^{14}\\
& + 2 e^2 m^6 n^{12} + 5 e^2 m^8 n^{10} + 2 e^2 m^{10} n^8 - e^2 m^{12} n^6 - e^4 m^2 n^{14} + 6 e^4 m^4 n^{12} + 15 e^4 m^6 n^{10}\\
& + 12 e^4 m^8 n^8 - 2 e^4 m^{10} n^6 - e^6 m^2 n^{12} + 9 e^6 m^4 n^{10} + 16 e^6 m^6 n^8 + 2 e^6 m^8 n^6 - 2 e^8 m^2 n^{10}\\
& + 9 e^8 m^4 n^8 + 4 e^8 m^6 n^6 + e^8 m^8 n^4 - 2 e^{10} m^2 n^8 + 2 e^{10} m^4 n^6 - m^6 n^{16} - 2 m^8 n^{14} + m^{12} n^{10}\\
& - 2 e^2 m^4 n^{16} - 4 e^2 m^6 n^{14} + 3 e^2 m^8 n^{12} + 2 e^2 m^{10} n^{10} - 2 e^4 m^4 n^{14} + 6 e^4 m^6 n^{12} + 7 e^4 m^8 n^{10}\\
& - e^4 m^{10} n^8 - 2 e^4 m^{12} n^6 - 2 e^6 m^2 n^{14} + 2 e^6 m^4 n^{12} + 13 e^6 m^6 n^{10} + 8 e^6 m^8 n^8 - 3 e^6 m^{10} n^6\\
& - 3 e^8 m^2 n^{12} + 3 e^8 m^4 n^{10} + 8 e^8 m^6 n^8 + 2 e^8 m^8 n^6 - 2 e^{10} m^2 n^{10} - 3 e^{10} m^4 n^8 - e^{12} m^2 n^8\\
& - m^8 n^{16} - m^{10} n^{14} - e^2 m^4 n^{18} - 3 e^2 m^6 n^{16} - 6 e^2 m^8 n^{14} - e^2 m^{10} n^{12} - 3 e^2 m^{12} n^{10}\\
& + e^4 m^4 n^{16} - 10 e^4 m^6 n^{14} - 5 e^4 m^8 n^{12} - 7 e^4 m^{10} n^{10} - 3 e^4 m^{12} n^8 - 7 e^6 m^4 n^{14} + 2 e^6 m^6 n^{12}\\
& - e^6 m^8 n^{10} - 2 e^6 m^{10} n^8 - 2 e^6 m^{12} n^6 - 2 e^8 m^2 n^{14} - 4 e^8 m^4 n^{12} + 2 e^8 m^6 n^{10} + 3 e^8 m^8 n^8\\
& - 2 e^8 m^{10} n^6 - 2 e^{10} m^2 n^{12} - 6 e^{10} m^4 n^{10} - 3 e^{10} m^6 n^8 - 3 e^{12} m^4 n^8 - 3 m^{10} n^{16} - m^{12} n^{14}\\
& - m^{14} n^{12} - e^2 m^6 n^{18} - 7 e^2 m^8 n^{16} - 10 e^2 m^{10} n^{14} - 5 e^2 m^{12} n^{12} - 2 e^2 m^{14} n^{10} - 15 e^4 m^6 n^{16}\\
& - 33 e^4 m^8 n^{14} - 22 e^4 m^{10} n^{12} - 13 e^4 m^{12} n^{10} - e^4 m^{14} n^8 - 2 e^6 m^4 n^{16} - 21 e^6 m^6 n^{14}\\
& - 21 e^6 m^8 n^{12} - 16 e^6 m^{10} n^{10} - 4 e^6 m^{12} n^8 - 11 e^8 m^4 n^{14} - 15 e^8 m^6 n^{12} - 15 e^8 m^8 n^{10}\\
& - 6 e^8 m^{10} n^8 - 2 e^8 m^{12} n^6 - e^{10} m^2 n^{14} - 9 e^{10} m^4 n^{12} - 8 e^{10} m^6 n^{10} - 5 e^{10} m^8 n^8 - e^{10} m^{10} n^6\\
& - 3 e^{12} m^6 n^8 - e^{14} m^4 n^8 + \cO([emn]^{28})\,.
\end{align*}

It is worth noting that the multi-graded HS lacks certain properties of the ungraded HS, such as a matching number of factors in the denominator with physical observables and a palindromic form in the numerator. In addition, as already mentioned in footnote <ref>, there is ambiguity when determining the form of the ungraded HS, similar ambiguity can arise for the multi-graded HS. For instance, it is possible to introduce a common factor $(1+e^2n^2)$ to both the numerator and denominator, then the factors $(1-e^2n^2)^2$ in denominator become $(1-e^2n^2)(1-e^4n^4)$, the numerator will also change accordingly. The form of the HS is less constrained in the multi-graded case. This ambiguity prevents us from further exploring the invariant structures of the theory. Instead of decoding the HS, analyzing the PL is more helpful in our case, since it is unique for both ungraded and multi-graded HS.

The invariants in our theory form a non-complete intersection ring making the PL a non-terminating series as a result. According to Eq. (<ref>), the PL can only be calculated up to some given order in the spurions. However, as we discussed in Sec. <ref>, we are only interested in the terms in the PL up to the first purely negative order, which we find to occur at order $\cO([emn]^{26})$. The PL up to this order is given by

\begin{align*}
\label{eq:gradedPL}
& \PL(e,m,n)= \( e^2+m^2+n^2\) + \( e^4+m^4+n^4+e^2 n^2+m^2 n^2\) + \( e^6+m^6+n^6+e^2 n^4 \right. \\
& \left. + e^4 n^2+2 m^2 n^4+m^4 n^2+e^2 m^2 n^2 \) + \( e^4 n^4+m^2 n^6+3 m^4 n^4+3 e^2 m^2 n^4 + e^2 m^4 n^2  \right. \\
& \left. + e^4 m^2 n^2 \) + \( m^2 n^8+2 m^4 n^6+2 m^6 n^4+5 e^2 m^2 n^6+4 e^2 m^4 n^4 + 4 e^4 m^2 n^4+e^4 m^4 n^2 \) \\
& + \( e^6 n^6+3 m^4 n^8+3 m^6 n^6+m^8 n^4+5 e^2 m^2 n^8 + 9 e^2 m^4 n^6+3 e^2 m^6 n^4 +7 e^4 m^2 n^6 \right. \\
& \left. +6 e^4 m^4 n^4+2 e^6 m^2 n^4\) + \( m^4 n^{10} + 3 m^6 n^8+m^8 n^6+3 e^2 m^2 n^{10}+11 e^2 m^4 n^8 \right. \\
& \left. +8 e^2 m^6 n^6+e^2 m^8 n^4+9 e^4 m^2 n^8+14 e^4 m^4 n^6+5 e^4 m^6 n^4+6 e^6 m^2 n^6+3 e^6 m^4 n^4 \right. \\
& \left. +e^8 m^2 n^4 \)+ \( m^4 n^{12}+m^6 n^{10}+e^2 m^2 n^{12}+8 e^2 m^4 n^{10}+8 e^2 m^6 n^8+3 e^2 m^8 n^6 \right. \\
& \left. +6 e^4 m^2 n^{10}+20 e^4 m^4 n^8+13 e^4 m^6 n^6+2 e^4 m^8 n^4+8 e^6 m^2 n^8 +13 e^6 m^4 n^6 \right. \\
& \left. +3 e^6 m^6 n^4+3 e^8 m^2 n^6+2 e^8 m^4 n^4 \) + \(2 e^2 m^4 n^{12}+3 e^4 m^2 n^{12}+11 e^4 m^4 n^{10} \right. \\
& \left. +11 e^4 m^6 n^8+5 e^4 m^8 n^6+5 e^6 m^2 n^{10} +18 e^6 m^4 n^8+13 e^6 m^6 n^6+e^6 m^8 n^4+4 e^8 m^2 n^8 \right. \\
& \left. +8 e^8 m^4 n^6+2 e^8 m^6 n^4 +e^{10} m^2 n^6-3 m^8 n^{10}-3 m^{10} n^8-3 e^2 m^6 n^{10}-4 e^2 m^8 n^8 \right. \\
& \left. -e^2 m^{10} n^6 \) + \( 2 e^6 m^6 n^8 + 4 e^6 m^8 n^6 + 9 e^8 m^4 n^8 + 7 e^8 m^6 n^6 + e^8 m^8 n^4 + 3 e^{10} m^4 n^6 \right. \\
& \left. - m^6 n^{14} - 9 m^8 n^{12} - 8 m^{10} n^{10} - 4 m^{12} n^8 - 2 e^2 m^4 n^{14} - 21 e^2 m^6 n^{12} - 32 e^2 m^8 n^{10} \right. \\
& \left. - 14 e^2 m^{10} n^8 - e^2 m^{12} n^6 - 9 e^4 m^4 n^{12} - 28 e^4 m^6 n^{10} - 18 e^4 m^8 n^8 - 2 e^4 m^{10} n^6 \right. \\
& \left. - e^6 m^4 n^{10} \) + \( 3 e^8 m^8 n^6 + 2 e^{10} m^6 n^6 - m^6 n^{16} - 10 m^8 n^{14} - 18 m^{10} n^{12} - 9 m^{12} n^{10} \right. \\
& \left. - 2 m^{14} n^8 - 2 e^2 m^4 n^{16} - 34 e^2 m^6 n^{14} - 76 e^2 m^8 n^{12} - 55 e^2 m^{10} n^{10} - 12 e^2 m^{12} n^8 \right. \\