1 Introduction

The observation in 2012 of a scalar particle with 125 GeV by the ATLAS and CMS collaborations [1, 2] has incentivized experimental searches for beyond the Standard Model (SM) particles at the LHC. On par with these experimental endeavors, theoretical efforts in the search for extra scalar particles have been strengthened since this discovery. A promising framework is found in N-Higgs doublet models (NHDM).

Such models have many free parameters, which are often curtailed by imposing some discrete family symmetry. Here, we focus on the implementation of \(A_4\) in a three Higgs doublet model (3HDM). The \(A_4\) group is the group of even permutations on 4 elements. It is the smallest discrete group to contain a three-dimensional irreducible representation (irrep), which is ideal for describing the three families of quarks with a minimal number of independent Yukawa couplings. Thus, NHDM supplemented by the \(A_4\) discrete symmetry has long been of interest in flavour physics research.

A number of early articles include: [3], mainly devoted to the leptonic sector and where the solution to the quark sector is briefly mentioned to include a fourth Higgs doublet and all quark fields in singlets (which is effectively the same as the Standard Model quark sector); [4], where \(A_4\) is broken by dimension four Yukawa couplings, which, upon renormalization, will affect the scalar potential [5], which requires three Higgs doublets in the down-type quark sector and a further two in the up-type quark sector, consisting of a 5HDM; and [6], which is devoted to the leptonic sector, but has the interesting side query that it might be possible to recover a realistic CKM matrix through soft-breaking of \(A_4\).

Quark mass matrices in the context of a 3HDM with Higgs doublets in the triplet representation of \(A_4\) were studied in [7, 8], with the vacuum expectation value (vev) structure \((e^{i\alpha }, e^{- i \alpha }, r)\), where \(\alpha \) and r are real constants. This vacuum solution was also included in the original study of the \(A_4\)-3HDM vacua in Ref. [9]. Unfortunately, Degee et al. [10] proved in 2013 that such a vacuum can never be the global minimum of the \(A_4\) symmetric 3HDM. In this beautiful paper, geometric techniques were used in order to identify all possible global minima (thus, all possible viable vacua) of the \(A_4\) symmetric 3HDM. Immediately thereafter, those minima were used to show that all assignments of the quark fields into irreps of \(A_4\), when combined with the possible vevs for the exact \(A_4\) potential, yield vanishing quark masses and/or a CP conserving CKM matrix, both of which are forbidden by experiment. This is in fact a consequence of a much broader theorem, proved in [11, 12]: given any flavour symmetry group, one can obtain a physical CKM mixing matrix and, simultaneously, non-degenerate and non-zero quark masses only if the vevs of the Higgs fields break completely the full flavour group. The idea is that a symmetry will reduce the number of redundant Yukawa couplings present in the SM, and it might even predict relations among observables which turn out to be consistent with experiment.

When studying in detail the extensions of \(A_4\) to the quark sector found by Ref. [13], we noticed that, in some of them, if it weren’t for the particular form of the vevs allowed by the exact \(A_4\) 3HDM potential, the Yukawa matrices could allow for massive quarks, and for a realistic CKM matrix. Since the \(A_4\) symmetric potential doesn’t allow for minima other than those shown in [10], here we consider the case where the \(A_4\) symmetry is softly broken by the addition of quadratic terms to the potential. Such terms do not spoil the theory’s renormalizability, but break the \(A_4\) symmetry.

Our article is organized as follows. We define the notation for the scalar potential in Sect. 2.1, discuss the Yukawa Lagrangian and the form of the possible mass matrices in Sect. 2.2, giving all the expressions needed for the fit in Sect. 2.3. In Sect. 3 we present our fit to the quarks mass matrices, while in Sect. 4 we discuss the viability of the vacuum found in the fit in terms of the scalar potential. Section 5 is devoted to the implementation of the theoretical constraints to be imposed, and in Sect. 6 we briefly discuss the constraints coming from the LHC. The results and conclusions are presented in Sects. 7 and 8, respectively. The Appendices contain some additional expressions that are needed for the fits.

2 Parameterization for the softly-broken \(A_4\) 3HDM

2.1 Potential and candidates for local minimum

The softly-broken potential of the 3HDM with an \(A_4\) symmetry is given by

$$\begin{aligned} V_H=V_{4,\, A_4}+M^2_{ij}\left( \phi _i^{\dagger } \phi _j \right) , \end{aligned}$$
(1)

where \(V_{4,\, A_4}\) is the quartic potential for the \(A_4\) symmetric three Higgs doublet model (3HDM), which is, in the notation of [10],

$$\begin{aligned} V_{4,\, A_4}&=\frac{\Lambda _0}{3} \left( \phi _1 ^ \dagger \phi _1 + \phi _2 ^ \dagger \phi _2 + \phi _3 ^ \dagger \phi _3\right) ^ 2\nonumber \\&\quad + \Lambda _1 \left[ \left( \text {Re}\left\{ \phi _1^\dagger \phi _2\right\} \right) ^2 + \left( \text {Re}\left\{ \phi _2^\dagger \phi _3\right\} \right) ^2\right. \nonumber \\&\quad \left. + \left( \text {Re}\left\{ \phi _3^\dagger \phi _1\right\} \right) ^2 \right] + \Lambda _2 \left[ \left( \text {Im}\left\{ \phi _1^\dagger \phi _2\right\} \right) ^2\right. \nonumber \\&\quad \left. + \left( \text {Im}\left\{ \phi _2^\dagger \phi _3\right\} \right) ^2 + \left( \text {Im}\left\{ \phi _3^\dagger \phi _1\right\} \right) ^2 \right] \nonumber \\&\quad + \frac{\Lambda _3}{3} \left[ (\phi _1 ^ \dagger \phi _1) ^ 2 + (\phi _2 ^ \dagger \phi _2) ^ 2 + (\phi _3 ^ \dagger \phi _3) ^ 2\right. \nonumber \\&\quad \left. - (\phi _1^\dagger \phi _1) (\phi _2^\dagger \phi _2) - (\phi _2^\dagger \phi _2) (\phi _3^\dagger \phi _3) - (\phi _3^\dagger \phi _3) (\phi _1^\dagger \phi _1)\right] \nonumber \\&\quad + \Lambda _4 \left[ \text {Re}\left\{ \phi _1^\dagger \phi _2\right\} \, \text {Im}\left\{ \phi _1^\dagger \phi _2\right\} + \text {Re}\left\{ \phi _2^\dagger \phi _3\right\} \,\right. \nonumber \\&\quad \times \left. \text {Im}\left\{ \phi _2^\dagger \phi _3\right\} + \text {Re}\left\{ \phi _3^\dagger \phi _1\right\} \, \text {Im}\left\{ \phi _3^\dagger \phi _1\right\} \right] . \end{aligned}$$
(2)

The matrix \(M^2_{ij}\) is a general hermitian matrix, which can be parameterized by

$$\begin{aligned} (M^2_{ij}) = \begin{pmatrix} m^2_{11} &{} m^2_{12} e^{i\theta _{12}} &{} m^2_{13} e^{i\theta _{13}}\\ m^2_{12} e^{-i\theta _{12}} &{} m^2_{22} &{} m^2_{23} e^{i\theta _{23}}\\ m^2_{13} e^{-i\theta _{13}} &{} m^2_{23} e^{-i\theta _{23}} &{} m^2_{33} \end{pmatrix}, \end{aligned}$$
(3)

where \(m^2_{ij}\) are real parameters with the dimension of mass squared.Footnote 1

Additionally, in the notation of [14], the exact \(A_4\) potential can be written as

$$\begin{aligned} V_{A_4}&= \frac{r_1 + 2 r_4}{3} \left[ (\phi _1^\dagger \phi _1) + (\phi _2^\dagger \phi _2) + (\phi _3^\dagger \phi _3) \right] ^2\nonumber \\&\quad + \frac{2 (r_1 - r_4)}{3} \left[ (\phi _1^\dagger \phi _1)^2 + (\phi _2^\dagger \phi _2)^2 \right. \nonumber \\&\quad \left. + (\phi _3^\dagger \phi _3)^2 - (\phi _1^\dagger \phi _1) (\phi _2^\dagger \phi _2) - (\phi _2^\dagger \phi _2) (\phi _3^\dagger \phi _3)\right. \nonumber \\&\quad \left. - (\phi _3^\dagger \phi _3) (\phi _1^\dagger \phi _1)\right] \nonumber \\&\quad + 2 r_7 \left( |\phi _1^\dagger \phi _2|^2 + |\phi _2^\dagger \phi _3|^2 + |\phi _3^\dagger \phi _1|^2 \right) \nonumber \\&\quad + \Big [ c_3 \left[ (\phi _1^\dagger \phi _2)^2 + (\phi _2^\dagger \phi _3)^2 + (\phi _3^\dagger \phi _1)^2 \right] + h.c. \Big ] . \end{aligned}$$
(4)

The relation between the two notations is

$$\begin{aligned}&r_1=\frac{1}{3}(\Lambda _0+\Lambda _3) ,\quad r_4=\frac{1}{6}(2 \Lambda _0-\Lambda _3) , \nonumber \\&r_7{=}\frac{1}{4}(\Lambda _1{+}\Lambda _2) , \nonumber \\&\text {Re}(c_3)= \frac{1}{4}(\Lambda _1-\Lambda _2) ,\quad \text {Im}(c_3)= -\frac{1}{4} \Lambda _4 . \end{aligned}$$
(5)

We consider that the scalar fields can take complex vacuum expectation values (vevs), to be determined later. Thus, we write,

$$\begin{aligned} \phi _i= \begin{bmatrix} \varphi ^+_i\\ \frac{|v_i| e^{i \rho _i}}{\sqrt{2}} + \frac{1}{\sqrt{2}} \left( x_i + i x_{i+3} \right) \end{bmatrix}. \end{aligned}$$
(6)

Because CP is spontaneously violated, the unrotated neutral fields have no definite CP, and for convenience we label them \(x_i,i=1,\ldots ,6\). We can also use the gauge freedom to absorb one of the phases in the vevs, that we choose to be \(\rho _1\). Therefore we have the vector of vevs defined as

$$\begin{aligned} \vec {v} = (|v_1|,|v_2| e^{i\rho _2},|v_3| e^{i\rho _3}). \end{aligned}$$
(7)

This vev contributes with four free parameters to our model, because one of the parameters is constrained by the mass of the gauge bosons to match the observed SM values,

$$\begin{aligned} |v_1|^2+ |v_2|^2+ |v_3|^2 \equiv v^2 \simeq (246\,\text {GeV})^2. \end{aligned}$$
(8)
Table 1 Extensions of \(A_4\) to the Yukawa sector with non-vanishing determinant, and non-zero J for general, complex valued, vevs \((v_1,v_2,v_3)\). In the table, \(\text {I}_d\) stands for the matrix \(M_d\) for case \(\text {I}\) and similarly for the other entries

The vev can also be parameterized as

$$\begin{aligned} \vec {v} = v \left( \cos (\beta _1)\cos (\beta _2), \cos (\beta _2)\sin (\beta _1) e^{ip_2}, \sin (\beta _2)e^{ip_3}\right) .\qquad \end{aligned}$$
(9)

Of the quantities arising out of the scalar potential, the vevs are the only relevant to the quark mass matrices. This leads many authors to just proclaim some vevs, without checking whether they can indeed be the global minima of a realistic Higgs potential. We will perform this crucial verification below, in Sect. 4.

2.2 Yukawa Lagrangian

As in Refs. [7, 13], we consider that the Higgs doublets are in the \({\textbf {3}}\) of \(A_4\) as well as the three left-handed SU(2) doublets \(Q_{Lj}\) of hypercharge 1/6. There are three right-handed SU(2) singlets \(n_{R,j}\) of hypercharge \(-1/3\) and three right-handed SU(2) singlets \(p_{R,j}\) of hypercharge 2/3. Our assignments for the singlets are as follows

$$\begin{aligned} n_{R1}, p_{R1}\rightarrow {\textbf {1}}, \quad n_{R2}, p_{R2} \rightarrow {\textbf {1}}^\prime , \quad n_{R3}, p_{R3} \rightarrow {\textbf {1}}^{\prime \prime } \ \text { of } A_4. \end{aligned}$$
(10)

Then, the \(A_4\) transformations on the fields are generated by [7, 13]

$$\begin{aligned} T: {\left\{ \begin{array}{ll} \phi _1\rightarrow \phi _2\rightarrow \phi _3\rightarrow \phi _1,\\ Q_{L1}\rightarrow Q_{L2}\rightarrow Q_{L3}\rightarrow Q_{L1},\\ n_{R1}\rightarrow n_{R1}, n_{R2}\rightarrow \omega n_{R2}, n_{R3}\rightarrow \omega ^2 n_{R3},\\ p_{R1}\rightarrow n_{R1}, p_{R2}\rightarrow \omega p_{R2}, p_{R3}\rightarrow \omega ^2 p_{R3}, \end{array}\right. } \end{aligned}$$
(11)

and

$$\begin{aligned} S: {\left\{ \begin{array}{ll} \phi _1\rightarrow \phi _1, \phi _2\rightarrow -\phi _2, \phi _3\rightarrow -\phi _3,\\ Q_{L1}\rightarrow Q_{L1}, Q_{L2}\rightarrow -Q_{L2}, Q_{L3}\rightarrow -Q_{L3}. \end{array}\right. } \end{aligned}$$
(12)

One can easily verify that the scalar potential in Eq. (4) is invariant under the previous transformations. Now we write the \(A_4\) invariant Yukawa Lagrangian for quarks. We have

$$\begin{aligned}&-{\mathcal {L}}_{\textrm{Yukawa}}= \sqrt{2}\, \hat{a}\left( \overline{Q_{L1}}\phi _1+ \overline{Q_{L2}}\phi _2+\overline{Q_{L3}}\phi _3 \right) n_{R1}\nonumber \\&\quad +\sqrt{2}\, \hat{b}\left( \overline{Q_{L1}}\phi _1+ \omega \, \overline{Q_{L2}}\phi _2 +\omega ^2\, \overline{Q_{L3}}\phi _3 \right) n_{R2}\nonumber \\&\quad +\sqrt{2}\, \hat{c}\left( \overline{Q_{L1}}\phi _1+ \omega ^2\,\overline{Q_{L2}}\phi _2 +\omega \, \overline{Q_{L3}}\phi _3 \right) n_{R3}\nonumber \\&\quad +\sqrt{2}\, \hat{a}'\left( \overline{Q_{L1}}\tilde{\phi }_1+ \overline{Q_{L2}}\tilde{\phi }_2+ \overline{Q_{L3}}\tilde{\phi }_3 \right) p_{R1}\nonumber \\&\quad +\sqrt{2}\, \hat{b}'\left( \overline{Q_{L1}}\tilde{\phi }_1+ \omega \, \overline{Q_{L2}}\tilde{\phi }_2 +\omega ^2\, \overline{Q_{L3}}\tilde{\phi }_3 \right) p_{R2}\nonumber \\&\quad +\sqrt{2}\, \hat{c}'\left( \overline{Q_{L1}}\tilde{\phi }_1+ \omega ^2\,\overline{Q_{L2}}\tilde{\phi }_2 +\omega \, \overline{Q_{L3}}\tilde{\phi }_3 \right) p_{R3} + \text {h.c.} , \end{aligned}$$
(13)

where, as usual,

$$\begin{aligned} \tilde{\phi }_j \equiv i\, \sigma _2 \phi ^*_j, \end{aligned}$$
(14)

and we define

$$\begin{aligned}{} & {} \hat{a}{=}a e^{i\, \alpha },\ \hat{b}{=}b e^{i\, \beta },\ \hat{c}{=}c e^{i\, \gamma },\ \hat{a}'{=}a' e^{i\, \alpha '},\ \hat{b}'{=}b' e^{i\, \beta '},\nonumber \\{} & {} \quad \hat{c}'{=}c' e^{i\, \gamma '}, \end{aligned}$$
(15)

where \(a,b,c,a',b',c'\) are real and positive. This choice of invariant Lagrangian corresponds to the case I identified in Ref. [13] (see the next section).

2.3 Yukawa matrices, masses and CKM

We aim to fit six quark masses and four CKM matrix elements to the currently accepted SM values for these observables. Therefore, we’re interested in softly-broken \(A_4\) symmetric models with up to ten parameters. Reference [13] has studied all of the possible extensions of \(A_4\) to the fermion sector. Using their results, we can check which of them can accommodate non-vanishing quark masses, CKM mixing angles and CP violation by considering a general vev \(\vec {v}\). We take the Jarlskog invariant as a measure of CP violation [15]. Out of all possibilities, we are left with five of them, which we list in Table 1. There, A are real constants, \(\Omega \) are constants in the \([0, 2\pi [\) interval, \(\omega =e^{i \frac{2 \pi }{3}}\) \((\omega ^3=1)\) and \(^T\) is the transpose of the matrix.

In the table above, we have used the convention where the quarks’ mass terms are written as

$$\begin{aligned} -{\mathcal {L}}_{\text {Yukawa}} \supset \overline{n_L} M_d n_R + \overline{p_L} M_u p_R + \text {h.c.}, \end{aligned}$$
(16)

where h.c. stands for the hermitian conjugate.

In the Yukawa sector, there are ten observables, six masses, three mixing angles and one Jarlskog invariant, therefore, we would prefer to look for a case with ten parameters, or less. All possible neutral vevs of the 3HDM are consistent with the parameterization in Eq. (9), which consists of four free parameters that we can fit; two angles, and two phases. Looking at the cases in Table 1, we will see that it is possible to reduce the number of free parameters by performing both basis transformations to right-handed quarks and global \(U(1)_Y\) rephasings, both of which have no effect on the physical predictions of the theory.

For case I, the down quark mass matrices read

$$\begin{aligned} M_d = \begin{pmatrix} a e^{i\alpha } v_1 &{} be^{i\beta }v_1 &{} ce^{i\gamma }v_1 \\ a e^{i\alpha } v_2 &{} \omega be^{i\beta }v_2 &{} \omega ^2 ce^{i\gamma }v_2 \\ a e^{i\alpha } v_3 &{} \omega ^2 be^{i\beta }v_3 &{} \omega ce^{i\gamma }v_3 \\ \end{pmatrix}\, = D_v W D_a D_{\alpha } , \end{aligned}$$
(17)

where (remember that the \(v_i\) are complex)

$$\begin{aligned}{} & {} D_{v}= \text {diag}(v_1, v_2, v_3), \,\, D_{a} = \text {diag}(a,b,c), \,\, \nonumber \\{} & {} D_{\alpha } = \text {diag}(e^{i \alpha },e^{i \beta },e^{i \gamma }), \,\, W = \begin{pmatrix} 1 &{} 1 &{} 1\\ 1 &{} \omega &{} \omega ^2 \\ 1 &{} \omega ^2 &{} \omega \end{pmatrix}. \end{aligned}$$
(18)

We see that we can perform a unitary transformation to the right-handed quarks that removes all three phases \(\alpha \), \(\beta \), \(\gamma \). The same holds for \(M_{u}\), by performing the substitution \(A \rightarrow A'\), \(\Omega \rightarrow \Omega '\) and \(v_i \rightarrow v_i^*\). We note that the case I matrices were also used by Ref. [16] as the mass matrices for the charged leptons.

In this work, we study this case, that corresponds to the Lagrangian in Eq. (13). Then, given that \(D_{\alpha }D_{\alpha }^{\dagger }={\mathbb {1}}\) and \(D_{a}D_{a}^{\dagger }=D_{a^2} = \text {diag}(a^2,b^2,c^2)\), we find

$$\begin{aligned} H_d\equiv & {} M_d M_d^{\dagger } = D_{v}S_d D_{v}^{\dagger },\nonumber \\ H_u\equiv & {} M_u M_u^{\dagger } = D_{v}^{\dagger }S_u D_{v}, \end{aligned}$$
(19)

where \(S_d=W D_{a^2} W^{\dagger }\) and \(a^2 \rightarrow a'^2\) for the up quark case. This matrix can now be explicitly written out using appropriate parameters as

$$\begin{aligned} S_d= \begin{pmatrix} \Sigma _d &{} Z_de^{i\phi _d} &{} Z_de^{-i\phi _d}\\ Z_de^{-i\phi _d} &{} \Sigma _d &{} Z_de^{i\phi _d}\\ Z_de^{i\phi _d} &{} Z_de^{-i\phi _d} &{} \Sigma _d \end{pmatrix}, \end{aligned}$$
(20)

where \(\Sigma _d\) and \(Z_d\) are real, and

$$\begin{aligned} \Sigma _d\equiv & {} a^2+b^2+c^2, \nonumber \\ Z_d\ e^{i\phi _d}\equiv & {} a^2+\omega ^2 b^2 + \omega c^2, \end{aligned}$$
(21)

with corresponding primes for the up case. For completeness, the specific forms for \(H_d\) and \(H_u\) found after using the parameterizations in Eqs. (9) and (21) are written in Appendix A. The eigenvalues of the matrices \(H_d\) and \(H_u\) will be fitted for the (square of the) quark masses, \((m_d^2,m_s^2,m_b^2)\) and \((m_u^2,m_c^2,m_t^2)\), respectively.

We now turn to the Cabibbo–Kobayashi–Maskawa (CKM) matrix. As found by Branco and Lavoura [17], the absolute values of the CKM matrix can be obtained through calculating the traces of appropriate powers of the matrices \(H_u\) and \(H_d\). They observe that

$$\begin{aligned} \text {Tr}\left( H_u^a H_d^b\right) \equiv L_{ab}=\sum _{k,i}U_{ki} (D_u^a)_{kk}(D_d^b)_{ii}, \end{aligned}$$
(22)

where \(U_{ki}=|V_{ki}|^2\) and V is the CKM matrix. The CKM matrix is unitary and therefore U only has four independent entries. Consequently, in order to compute U, it is only necessary to resort to

$$\begin{aligned} L_{11}&= U_{ki} (D_u)_{kk}(D_d)_{ii} , \nonumber \\ L_{12}&= U_{ki} (D_u)_{kk}(D_d^2)_{ii}, \nonumber \\ L_{21}&= U_{ki} (D_u^2)_{kk}(D_d)_{ii}, \nonumber \\ L_{22}&= U_{ki} (D_u^2)_{kk}(D_d^2)_{ii}. \end{aligned}$$
(23)

These equations are linear in \(U_{ik}\) and are, therefore, invertible for this variable. Thus, by picking \(U_{11}\), \(U_{21}\), \(U_{13}\), and \(U_{23}\) (respectively, \(U_{ud}\), \(U_{cd}\), \(U_{ub}\), and \(U_{cb}\)), we are able to obtain a unique solution for the magnitudes of the CKM elements as a function of \(L_{ab}\) and the quark masses. Namely,

$$\begin{aligned} U_{11}= & {} \left( {m_b}^2-{m_s}^2\right) \left( {m_c}^2-{m_t}^2\right) \frac{ a_{11}}{\text {det}},\nonumber \\ U_{21}= & {} \left( {m_b}^2-{m_s}^2\right) \left( {m_u}^2-{m_t}^2\right) \frac{ a_{21}}{\text {det}}, \nonumber \\ U_{13}= & {} \left( {m_d}^2-{m_s}^2\right) \left( {m_c}^2-{m_t}^2\right) \frac{ a_{13}}{\text {det}}, \nonumber \\ U_{23}= & {} \left( {m_d}^2-{m_s}^2\right) \left( {m_u}^2-{m_t}^2\right) \frac{ a_{23}}{\text {det}}, \end{aligned}$$
(24)

where

$$\begin{aligned}{} & {} a_{11} = {L_{11}} \left( {m_b}^2+{m_s}^2\right) \left( {m_c}^2+{m_t}^2\right) \nonumber \\{} & {} \quad {-}{L_{12}} \left( {m_c}^2{+}{m_t}^2\right) {-}{L_{21}} \left( {m_b}^2{+}{m_s}^2\right) {+}{L_{22}} \nonumber \\{} & {} \quad {+}m_b^2 \left( {-}m_c^2 m_t^2 \left( m_d^2{+}m_s^2\right) {-}m_s^2 m_u^2 \left( m_c^2{+}m_t^2\right) {+}m_s^2 m_u^4\right) \nonumber \\{} & {} \quad {+}m_c^2 m_d^2 m_t^2 \left( m_d^2{-}m_s^2\right) , \end{aligned}$$
(25)
$$\begin{aligned}{} & {} a_{21} {=} {-}{L_{11}}\left( {m_b}^2{+}{m_s}^2\right) \left( {m_t}^2{+}{m_u}^2\right) {+}{L_{12}}\left( {m_u}^2{+}{m_t}^2\right) \nonumber \\{} & {} \quad {+}{L_{21}}\left( {m_b}^2{+}{m_s}^2\right) {-}{L_{22}} \nonumber \\{} & {} \quad {+}m_b^2 \left( m_c^2 m_s^2 \left( m_t^2+m_u^2{-}m_c^2\right) {+}m_t^2 m_u^2 \left( m_d^2{+}m_s^2\right) \right) \nonumber \\{} & {} \quad {+}m_d^2 m_t^2 m_u^2 \left( m_s^2{-}m_d^2\right) , \end{aligned}$$
(26)
$$\begin{aligned}{} & {} a_{13} {=} {-}{L_{11}}\left( {m_d}^2{+}{m_s}^2\right) \left( {m_t}^2{+}{m_c}^2\right) {+}{L_{12}}\left( {m_c}^2{+}{m_t}^2\right) \nonumber \\{} & {} \quad {+}{L_{21}}\left( {m_d}^2{+}{m_s}^2\right) {-}{L_{22}} \nonumber \\{} & {} \quad {+}m_b^2 m_c^2 m_t^2 \left( m_d^2{+}m_s^2{-}m_b^2\right) {+}m_d^2 m_s^2 \left( m_c^2 \left( m_t^2{+}m_u^2\right) \right. \nonumber \\{} & {} \quad \left. {+}m_u^2 \left( m_t^2{-}m_u^2\right) \right) , \end{aligned}$$
(27)
$$\begin{aligned}{} & {} a_{23} {=} {L_{11}}\left( {m_d}^2{+}{m_s}^2\right) \left( {m_t}^2{+}{m_u}^2\right) {-}{L_{12}}\left( {m_u}^2{+}{m_t}^2\right) \nonumber \\{} & {} \quad {-}{L_{21}}\left( {m_d}^2{+}{m_s}^2\right) {+}{L_{22}} \nonumber \\{} & {} \quad + m_t^2 m_u^2 \left( m_b^4-m_b^2 \left( m_d^2+m_s^2\right) -m_d^2 m_s^2\right) \nonumber \\{} & {} \quad +m_c^4 m_d^2 m_s^2-m_c^2 m_d^2 m_s^2 \left( m_t^2+m_u^2\right) , \end{aligned}$$
(28)

and

$$\begin{aligned}{} & {} \text {det} = \left( {m_b}^2-{m_d}^2\right) \left( {m_c}^2-{m_u} ^2\right) \left( {m_d}^2-{m_s}^2\right) \nonumber \\{} & {} \quad \times \left( {m_u} ^2-{m_t}^2\right) \left( {m_b}^2-{m_s}^2\right) \left( {m_c}^2-{m_t}^2\right) . \end{aligned}$$
(29)

In these equations, the \(L_{ij}\) are obtained by evaluating the left hand side of Eq. (22). Finally, we note that knowing these four CKM magnitudes, we can determine the Jarlskog invariant [15], up to its sign. Thus, given some phase convention, we are also able to determine the phases of all CKM matrix elements.

3 The fit to the quark mass matrices

3.1 Parameters and observables

We would like to fit 10 observables (6 quark masses and 4 CKM parameters) with the 10 free parameters that we have in this model,

$$\begin{aligned} \beta _1, \beta _2, \rho _2, \rho _3, \Sigma _d, \Sigma _u, Z_d, Z_u, \phi _d, \phi _u. \end{aligned}$$
(30)

Notice that this is a huge improvement over the SM, where there are 18 complex Yukawa parameters. Similarly, in Ref. [4], there are 18 Yukawa couplings; in their notation \(h_1^{u,d}\), \(h_2^{u,d}\), \(h_3^{u,d}\), and those with \(h \rightarrow h'\) and \(h \rightarrow h''\). These reduce to 12 complex parameters, even after the approximation in their equation (19). So, having only 10 real parameters is already excellent.

Moreover, our 10 parameters are constrained. Although we were not able to find an analytical relation which expresses such a constraint, we can show numerically that it does exist. We postpone this proof until the end of Sect. 3.3. The upshot is that it was not guaranteed a priori that our 10 parameters would be able to fit the 10 observables. Turning the argument around, the fact that the 10 experimental values do allow for a good fit in the \(A_4\)-3HDM can be viewed as a success for the model.

3.2 The fitting procedure

We have implemented a \(\chi ^2\) analysis of the model, through a minimization performed using the CERN Minuit library [18]. The observables employed in this analysis, labeled by \(i=1,\ldots ,11\) are specified in Table 2, where \(\overline{X}_i\) represents the experimental mean value of the observable \(X_i\) and \(\sigma _i\) is the experimental error, which, when both left and right bounds are stated, is assumed to be the largest of the two.

Table 2 Experimental values and fit results

The data on the quark masses as well as for the CKM matrix elements and the Jarlskog invariant experimental values were obtained from [19]. As mentioned, |J| is fixed by \(|V_{11}|\), \(|V_{21}|\), \(|V_{13}|\), and \(|V_{23}|\). However, using it in the fit speeds the numerical convergence onto a good solution.

The \(\chi ^2\) function depends on the 10 parameters of our model (30),

$$\begin{aligned} \beta _1, \beta _2, \rho _2, \rho _3, \Sigma _d, \Sigma _u, Z_d, Z_u, \phi _d, \phi _u \end{aligned}$$
(31)

and is written as

$$\begin{aligned} \chi ^2({{\textbf {p}}})= \sum _{i=1}^{11} \left( \frac{P_i({{\textbf {p}}})-\overline{X}_i}{\sigma _i} \right) ^2, \end{aligned}$$
(32)

where \(P_i({{\textbf {p}}})\) is our model’s prediction for each of the 11 (10 + J) observables. The fit is complicated by the fact that the masses (squared) are obtained from the eigenvalues of \(H_d,H_u\) but the elements of the CKM also depend on the masses, see Eq. (24). So, we start by calculating the eigenvalues of \(H_d\) and \(H_u\), which depend only on the parameters in Eq. (31). Then, we evaluate the \(L_{ij}\) from the left hand side of Eq. (22), and finally the CKM elements are obtained from Eq. (24). In Appendix A we give the explicit expressions for the matrices \(H_d\) and \(H_u\).

3.3 Results of the fit

We have found an excellent fit of our model to the data, given in the second column of Table 2. This fit results in \(\chi ^2= 0.058\), for the parameters

$$\begin{aligned} \beta _1&=1.4260868\ \text {radians} ,\nonumber \\ \beta _2&= 1.5424328\ \text {radians} ,\nonumber \\ \rho _2&=4.2784971\ \text {radians} ,\nonumber \\ \rho _3&=5.3682785\ \text {radians} ,\nonumber \\ \Sigma _d&= 0.2889178\times 10^{-3} ,\nonumber \\ \Sigma _u&=0.4927455 ,\nonumber \\ Z_d&=0.1816577\times 10^{-3} ,\nonumber \\ Z_u&=0.4758317 ,\nonumber \\ \phi _d&=-1.7324779 \ \text {radians} ,\nonumber \\ \phi _u&=0.20644967\times 10^{-1} \ \text {radians} . \end{aligned}$$
(33)

This fit also leads to the data in the third column of Table 2, as well as to the vevs

$$\begin{aligned} |v_i| = \left( 1.00604, 6.90357, 245.901\right) \,(\text {GeV}). \end{aligned}$$
(34)

We notice that the vevs obey \(v_1<v_2\ll v_3\). This hierarchy of vevs is related to the hierarchy of the quark masses. This was also obtained in Ref. [7], although their model is not consistent, as their vev structure is not that of [10] for the symmetric \(A_4\) potential they consider.

We can now perform a second (toy) fitting procedure, which illustrates the fact that the ten parameters in our model are constrained, as announced at the end of Sect. 3.1. In this fit, we take all experimental values in Table 2, except that we trade the correct experimental value of \(m_s\) for \( m_s=(2 \pm 0.02)\,\text {GeV}\). Now, the fit is very poor, having \(\chi ^2=600\). If these had been the correct experimental values for the 10 observables, then our model would not be able to fit them. Conversely, the fact that such a fit is possible is a success for the model.

4 Viability of the vacuum found in the fit

We start by defining the three doublets as in Eq. (6). Next we define the physical eigenstates for the charged Higgs as \((G^+,S^+_2,S^2_3)^T\), and for the neutral states we have \((G^0,S^0_2,S^0_3,S^0_4,S^0_5,S^0_6)^T\), identifying the would-be Goldstone bosons \(G^+\equiv S^+_1\) and \(G^0\equiv S^0_1\). With these conventions, and following the definitions in [20], we define the \(3\times 3\) matrix \(\tilde{U}\) by

$$\begin{aligned} \varphi _i^+ \equiv \sum _{j=1}^3 \tilde{U}_{ij} S^+_j, \end{aligned}$$
(35)

and the \(3\times 6\) matrix \(\tilde{V}\) by

$$\begin{aligned} x_i + i x_{i+3} = \sum _{j=1}^6 \tilde{V}_{ij} S^0_j. \end{aligned}$$
(36)

These matricesFootnote 2 are then related to the diagonalization matrices of the charged and neutral scalars, to which we now turn.

4.1 The minimization of the potential

In our procedure we already know the values of the vevs. So, we use the stationarity equations to solve for the soft parameters, and leave the quartic parameters of the potential \(\Lambda _i\) as free parameters. In this way we can solve for \(m^2_{11}, m^2_{22}, m^2_{33}\) as well as for \(\text {Im} (m^2_{12}), \text {Im} (m^2_{13})\), leaving as free parameters the \(\Lambda _i\) and \(\text {Re} (m^2_{12}), \text {Re} (m^2_{13}), \text {Re} (m^2_{23}),\text {Im} (m^2_{23})\). When evaluating the scalar mass matrices (see below) the conditions have to be applied to ensure that we are at the minimum. For completeness we write these conditions in Appendix B.

4.2 The charged mass matrix

The charged mass matrix is obtained from the second derivatives at the minimum,

$$\begin{aligned} {\mathcal {M}}^2_C = \left. \frac{\partial ^2 V_H}{\partial \varphi _i^+\partial \varphi _j^-}\right| _{\text {Min}}. \end{aligned}$$
(37)

The matrix \({\mathcal {M}}^2_C\) is an hermitian matrix, with real eigenvalues and satisfying, with our usual conventions,

$$\begin{aligned} R_{\textrm{ch}} {\mathcal {M}}^2_C R^\dagger _{\textrm{ch}} =\text {diag}(0,m^2_{S^+_2},m^2_{S^+_3})\equiv {\mathcal {M}}^2_{D_{ch}}, \end{aligned}$$
(38)

where \(R_{\textrm{ch}}\) is an unitary matrix that satisfies,

$$\begin{aligned} S^+_i = \sum _{j=1}^3 \left( R_{\textrm{ch}}\right) _{ij} \varphi ^+_j. \end{aligned}$$
(39)

This can be seen from

$$\begin{aligned} {\mathcal {L}}_{\textrm{mass}}&= - \varphi ^{-}_i \left( {\mathcal {M}}^2_C\right) _{ij}\varphi ^{+}_j {=} {-} \varphi ^{-}_i \left( R^\dagger _{\textrm{ch}} R_{\textrm{ch}} {\mathcal {M}}^2_C R^\dagger _{\textrm{ch}} R_{\textrm{ch}} \right) _{ij}\varphi ^{+}_j \nonumber \\&= - \varphi ^{-}_i \left( R^\dagger _{\textrm{ch}} {\mathcal {M}}^2_{D_{ch}} R_{\textrm{ch}} \right) _{ij} \varphi ^{+}_j\nonumber \\&= - S^{-}_i \left( {\mathcal {M}}^2_{D_{ch}}\right) _{ij} S^+_j, \end{aligned}$$
(40)

where we have used Eq. (39).

We have checked both algebraically and numerically that we have a zero eigenvalue corresponding to \(G^+\) and we require that all other masses squared are positive, a condition for a local minimum.

4.3 The neutral mass matrix

Since in our case CP is not conserved, we denote the unrotated neutral scalars by \(x_i, i=1,\ldots ,6\), as in Eq. (6). We therefore obtain the neutral mass matrix as,

$$\begin{aligned} {\mathcal {M}}^2_N = \left. \frac{\partial ^2 V_H}{\partial x_i\partial x_j}\right| _{\text {Min}}. \end{aligned}$$
(41)

This is a symmetric real matrix diagonalized by an orthogonal \(6\times 6\) matrix,

$$\begin{aligned} R_{\textrm{neu}} {\mathcal {M}}^2_N R^T_{\textrm{neu}}{} & {} =\text {diag}(0,m^2_{S^0_2},m^2_{S^0_3},m^2_{S^0_4},m^2_{S^0_5},m^2_{S^0_6})\nonumber \\{} & {} \equiv {\mathcal {M}}^2_{D_{\textrm{neu}}}, \end{aligned}$$
(42)

with

$$\begin{aligned} S^0_i = \sum _{j=1}^6 \left( R_{\textrm{neu}}\right) _{ij} x_j. \end{aligned}$$
(43)

As for the case of the charged scalars, we have checked both algebraically and numerically that we have a zero eigenvalue corresponding to \(G^0\) and we require that all other masses squared are positive, a condition for a local minimum.

5 Theoretical constraints

After having shown that a solution exists for the vevs and parameters in the Yukawa sector that correctly fits the quarks masses and the CKM entries, we have to show that this is compatible with the scalar potential analysis. In particular we have to show that the vevs correspond to a local minimum of the potential and that both the theoretical constraints as well as those coming from LHC are satisfied. In this section we analyze the theoretical constraints.

5.1 Perturbative unitarity

This problem was already solved in [14], so we take the potential in the form of Eq. (4). From Ref. [14] we have the following expression for the eigenvalues \(\lambda _i\)Footnote 3

$$\begin{aligned} \lambda _1&=2 \left( 2 \text {Re}(c_3)+r_1\right) \end{aligned}$$
(44)
$$\begin{aligned} \lambda _2&=2 \left( \sqrt{3}\, |\text {Im}(c_3)| -\text {Re}(c_3) + r_1\right) \end{aligned}$$
(45)
$$\begin{aligned} \lambda _3&=2 \left( -\sqrt{3}\, |\text {Im}(c_3)| -\text {Re}(c_3) + r_1\right) \end{aligned}$$
(46)
$$\begin{aligned} \lambda _4&=2 (r_4+r_7) \end{aligned}$$
(47)
$$\begin{aligned} \lambda _5&=2 (r_4-r_7) \end{aligned}$$
(48)
$$\begin{aligned} \lambda _6&=2 (r_1+2 r_7) \end{aligned}$$
(49)
$$\begin{aligned} \lambda _7&=2 (r_1-r_7) \end{aligned}$$
(50)
$$\begin{aligned} \lambda _8&=2 (r_4+|c_3|) \end{aligned}$$
(51)
$$\begin{aligned} \lambda _9&=2 (r_4-|c_3|) \end{aligned}$$
(52)
$$\begin{aligned} \lambda _{10}&=6 r_1+8 r_4+4 r_7 \end{aligned}$$
(53)
$$\begin{aligned} \lambda _{11}&=6 r_1-2 (2 r_4+r_7) \end{aligned}$$
(54)
$$\begin{aligned} \lambda _{12}&=6 |c_3|+2 r_4+4 r_7 \end{aligned}$$
(55)
$$\begin{aligned} \lambda _{13}&=-6 |c_3|+2 r_4+4 r_7. \end{aligned}$$
(56)

Perturbative unitarity is satisfied if

$$\begin{aligned} |\lambda _i| < 8 \pi , \quad \forall i. \end{aligned}$$
(57)

5.2 The BFB conditions

For the \(A_4\) symmetric potential, the conditions for boundedness from below along the neutral directions (BFB-n) have been conjectured in [21], and proved to hold in [22]. These are

$$\begin{aligned}&\Lambda _0 + \Lambda _3 \ge 0, \end{aligned}$$
(58)
$$\begin{aligned}&\frac{4}{3}(\Lambda _0 + \Lambda _3) + \frac{1}{2} (\Lambda _1 + \Lambda _2) - \Lambda _3\nonumber \\&\quad - \frac{1}{2} \sqrt{(\Lambda _1 - \Lambda _2)^2+\Lambda _4^2} \ge 0, \end{aligned}$$
(59)
$$\begin{aligned}&\Lambda _0 + \frac{1}{2} (\Lambda _1 + \Lambda _2) + \frac{1}{2} (\Lambda _1 - \Lambda _2) \cos {(2 k \pi /3 )}\nonumber \\&\quad + \frac{1}{2} \Lambda _4 \sin {(2 k \pi /3 )} \ge 0\quad (k=1,2,3). \end{aligned}$$
(60)

However, as shown in [21, 23], a potential which is BFB-n is not necessarily BFB along the charge breaking directions (BFB-c). Necessary BFB-c conditions have yet to be found for the \(A_4\) 3HDM, but sufficient conditions have been proposed in [24] following the technique developed in [25]. They are,

$$\begin{aligned} A_d \ge 0,\quad A_o \ge -A_d/2, \end{aligned}$$
(61)

where

$$\begin{aligned} A_d&= a=\frac{2}{3}(\Lambda _0+\Lambda _3),\nonumber \\ A_o&= b + \text {min}(0,c)-d\nonumber \\&=\frac{1}{3}(2\Lambda _0 - \Lambda _3) +\frac{1}{2}(\Lambda _1+\Lambda _2) \nonumber \\&\quad +\min (0,-\frac{1}{2}(\Lambda _1+\Lambda _2)) - \frac{1}{2} \sqrt{(\Lambda _1-\Lambda _2)^2+\Lambda _4^2}. \end{aligned}$$
(62)

It is important to remark that, since these are sufficient, but not necessary, conditions, some good points in parameter space may be excluded by this restriction.

5.3 The oblique parameters STU

For this we use the notation and results from [20], which require the matrices \(\tilde{U}\) and \(\tilde{V}\). Comparing Eq. (39) with the definition in Eq. (35), we conclude that

$$\begin{aligned} \tilde{U} = R_{\textrm{ch}}^\dagger , \end{aligned}$$
(63)

where the matrix \(R_{\textrm{ch}}\) is obtained from the numerical diagonalization of Eq. (38). Similarly, comparing Eq. (43) with the definition of \(\tilde{V}\) in Eq. (36), we get,

$$\begin{aligned} \tilde{V}_{ij}= \left( R^T_{\textrm{neu}}\right) _{ij}+i\left( R^T_\textrm{neu}\right) _{i+3,j}. \end{aligned}$$
(64)

Having \(\tilde{U}\) and \(\tilde{V}\), we can construct the needed matrices \(\text {Im}\left( \tilde{V}^\dagger \tilde{V} \right) \), \(\tilde{U}^\dagger \tilde{U}\), \(\tilde{V}^\dagger \tilde{V}\) and \(\tilde{U}^\dagger \tilde{V}\), and implement the procedure of [20].

5.4 Global minimum

After finding a set of \(m_{i,j}\) and \(\Lambda _i\) which reproduce the vevs in Eq. (34) necessary for a good fit of the quark mass matrices, and after performing the previous theoretical checks on the scalar potential, we must still ensure that our minimum is indeed the global minimum. This step is almost never taken in studies of quark mass matrices, since there are no exact analytical formulae for it. Moreover, one must check that there are no lower minima both along the neutral directions and along the charge breaking directions. We follow the strategy discussed in Ref. [24]. Take a specific set of \(m^2_{ij}\) and \(\Lambda _i\). Then we parameterize the scalar doublets as [23, 24],

$$\begin{aligned}{} & {} \langle \phi _1 \rangle =\sqrt{r_1}\, \begin{pmatrix} 0\\ 1 \end{pmatrix}, \quad \langle \phi _2 \rangle =\sqrt{r_2}\, \begin{pmatrix} \sin (\alpha _2)\\ \cos (\alpha _2) e^{i \beta _2} \end{pmatrix}, \quad \nonumber \\{} & {} \langle \phi _3 \rangle =\sqrt{r_3} e^{i \gamma }\, \begin{pmatrix} \sin (\alpha _3)\\ \cos (\alpha _3) e^{i \beta _3} \end{pmatrix}, \end{aligned}$$
(65)

where we have already used the gauge freedom. Now we let the vevs run free, for both charge conserving and charge violating directions. We give one seed point and perform a minimization of the potential using the CERN Minuit library [18]. We obtain not only the value of the potential at the minimum, but also the values of \(r_i,\alpha _2,\beta _2,\alpha _3,\beta _3\) and \(\gamma \). Then, we take one more (randomly generated) seed point and repeat the minimization. Finally, we take the minimum as the global one if it is found as the global minimum in each of 200 searches with randomly generated seed points. We have done this verification for every point that passed all the constraints. In all cases, we found that the local minimum was also a global minimum. In particular we always found that

$$\begin{aligned} \sin (\alpha _2)=\sin (\alpha _3)=0, \end{aligned}$$
(66)

showing that we do not have charged breaking directionsFootnote 4 and, comparing with Eq. (6), we verified numerically that,

$$\begin{aligned}{} & {} \frac{|v_i|}{\sqrt{2}}= \sqrt{r_i},\quad e^{i\, \rho _2}= \cos (\alpha _2)\, e^{i\, \beta _2},\quad \nonumber \\{} & {} e^{i\, \rho _3}= \cos (\alpha _3)\, e^{i\, (\beta _3+\gamma )}. \end{aligned}$$
(67)

6 Simple LHC constraints

Up to now we have implemented the theoretical constraints on the model. The next step is to implement the LHC constraints. To do this completely one would have to implement all the decays of the neutral and charged Higgs as well as their branching ratios. One would also have to worry about the electric dipole moments (EDM) and the flavour-changing neutral couplings (FCNC), as the model does not have a structure of couplings of the Higgs to the fermions that automatically ensures vanishing FCNC [26,27,28]. This lies beyond the scope of the present work. Nonetheless, we can implement easily the constraints that come from \(h \rightarrow WW/ZZ\) in the \(\kappa \) formalism, where the deviation from the coupling of the SM Higgs boson to a pair of W’s (or Z’s) is measured by \(\kappa _V\).

Fig. 1
figure 1

Left panel: Relation between \(m_{S_3^0}\) and \(m_{S_4^0}\). Right panel: Relation between \(m_{S_3^0}\) and \(m_{H_1^+}\). Color conventions: No cuts (red); with cuts (green)

Fig. 2
figure 2

Left panel: Relation between \(m_{S_5^0}\) and \(m_{S_6^0}\). Right panel: Relation between \(m_{S_5^0}\) and \(m_{H_2^+}\). Color conventions: No cuts (red); with cuts (green)

Fig. 3
figure 3

Left panel: Relation between \(m_{S_3^0}\) and \(m_{S_5^0}\). Right panel: Relation between \(m_{H_1^+}\) and \(m_{H_2^+}\). Color conventions: No cuts (red); with cuts (green)

Fig. 4
figure 4

Left panel: Relation between \(\kappa _V\) and \(\Lambda _1\). Right panel: Relation between \(\kappa _V\) and \(\Lambda _4\). Color conventions: No cuts (red), with theoretical cuts (green), and after the \(\kappa _V\) constraint (blue)

Fig. 5
figure 5

Left panel: Relation between \(\Lambda _0\) and \(\Lambda _4\). Right panel: Relation between \(\Lambda _0\) and \(\Lambda _3\) Color conventions: No cuts (red), with theoretical cuts (green), and after the \(\kappa _V\) constraint (blue)

In our model,

$$\begin{aligned} \kappa _V= & {} R^{\textrm{neu}}_{21} v_1 + R^{\textrm{neu}}_{22} v_2 \cos (\rho _2) + R^{\textrm{neu}}_{23} v_3 \cos (\rho _3)\nonumber \\{} & {} + R^{\textrm{neu}}_{25} v_2 \sin (\rho _2) + R^{\textrm{neu}}_{26} v_3 \sin (\rho _3), \end{aligned}$$
(68)

where \(R^{\textrm{neu}}\) is matrix defined in Eq. (42). We take the experimental constraint from ATLAS [29],

$$\begin{aligned} \kappa _W= 1.0206{\,}^{+ 0.05172}_{- 0.05087},\quad \kappa _Z= 0.99{\,}^{+ 0.06136}_{- 0.05214}. \end{aligned}$$
(69)

7 Results

In this section we present the results of the analysis of the scalar potential after imposing that we have a good solution for the fit of the quarks masses and CKM entries, as explained in Sect. 3.

7.1 Scanning strategy

We start by imposing the vevs obtained in the fit.

$$\begin{aligned}{} & {} v_1=1.00604\, (\text {GeV}),\quad v_2=6.90357\, e^{i\, 4.278497}\, (\text {GeV}), \quad \nonumber \\{} & {} v_3=245.901\, e^{i\, 5.368278} \, (\text {GeV}). \end{aligned}$$
(70)

Now we vary the free parameters of the potential in the following ranges,

$$\begin{aligned}{} & {} \log _{10}|\Lambda _i| \in [-3,1], \quad \log _{10}|\text {Im}(m^2_{23})| \in [-1,7]\,\text {GeV}^2, \quad \nonumber \\{} & {} \log _{10}|\text {Re}(m^2_{ij})| \in [-1,7]\,\text {GeV}^2, \end{aligned}$$
(71)

where in the last equation we use

$$\begin{aligned} m^2_{ij} \in \left\{ m^2_{12},m^2_{13},m^2_{23} \right\} . \end{aligned}$$
(72)

We randomly scan as in Eq. (71), and then:

  1. 1.

    Apply the theoretical constraints that only depend on the \(\Lambda _i\), that is BFB and perturbative unitarity.

  2. 2.

    Then obtain the eigenvalues for the charged and neutral scalars. Verify that all the masses squared are positive, and that we have a zero eigenvalue corresponding to the Goldstone bosons, \(G^0\) and \(G^+\).

  3. 3.

    Verify the S, T and U oblique parameters.

  4. 4.

    Apply the LHC constraint on \(\kappa _V\).

  5. 5.

    Check numerically that the vev is indeed a global minimum.

7.2 The scalar spectrum

We found that there is a strong correlation in the scalar masses. If we denote the masses of the neutral scalars by \((m_{G^0}=0,m_{S^0_{2}},m_{S^0_{3}},m_{S^0_{4}},m_{S^0_{5}},m_{S^0_{6}})\), and \((m_{G^+}=0,m_{H^+_{1}},m_{H^+_{2}})\) for the charged scalars, we find numerically that

$$\begin{aligned} m_{S^0_{3}}\simeq m_{S^0_{4}}\simeq m_{H^+_{1}},\quad m_{S^0_{5}}\simeq m_{S^0_{6}}\simeq m_{H^+_{2}}. \end{aligned}$$
(73)

This is true even if we do not require \(m_{S^0_{2}}=125\) GeV, and specially true after implementing the constraints of perturbative unitarity, BFB and STU. But, as we want to reproduce the LHC results, we also required that [19]

$$\begin{aligned} m_{S_2^0}= 125.25 \pm 0.17\,\text {GeV}. \end{aligned}$$
(74)

In the following figures we show the correlation among the masses. Included in red are the points generated before the theoretical cuts were applied, and in green the points remaining after the constraints were implemented (Figs. 1, 2, 3).

7.3 The \(\kappa _V\) constraint

We can now implement the \(\kappa _V\) constraint on the model. In the following figures, in red are points without cuts, in green with cuts but no \(\kappa _V\) constraint, and finally in blue points remaining after this constraint is applied. We took the ATLAS result of Eq. (69) at \(2\sigma \). While the theoretical constraints cut around 88% of the points, the \(\kappa _V\) constraint only cuts 22% of the remaining points. In Fig. 4 we show the relation between \(\kappa _V\) and \(\Lambda _{1,4}\) for the three sets of points as discussed above.

In fact it is not obvious from Fig. 4 that the \(\kappa _V\) constraint only cuts about 22% of the points that pass the other cuts. This is because there is a very large number of points with \(|\kappa _V|\lesssim 1\), even without theoretical cuts, and this is even more so after imposing the theoretical cuts. In this figure, we have 200,000 points in the green region, but from these 156,516 are in the blue region. That is, after theoretical cuts, 78% of the points also satisfy the \(\kappa _V\) constraint. In Fig. 5 we show the relation between \(\Lambda _0\) and \(\Lambda _{3,4}\) for the same sets of points. We see that, while for \((\Lambda _0,\, \Lambda _4)\) there is not much difference before and after the \(\kappa _V\) constraint, the same is not true for \((\Lambda _0,\, \Lambda _3)\), where the constraints impose a linear relation between those two parameters. We note that, while \(\Lambda _0\) is always positive, \(\Lambda _3\) can be negative respecting the BFB condition in Eq. (58), \(\Lambda _0+\Lambda _3 \ge 0\), as it is clear in the right-handed panel of Fig. 5. Before we end this section, let us remark that we did not redraw the figures in Sect. 7.2 after imposing the \(\kappa _V\) constraint, as the blue points would just superimpose the green points, as we have checked.

8 Conclusions

It is known that the 3HDM symmetric under an exact \(A_4\) symmetry is not compatible with non-zero quark masses and/or non-block-diagonal CKM matrix [13]. In this work, we studied a 3HDM with \(A_4\) softly broken. This allows us to evade the above result, by enlarging the structure of the possible vacua.

We obtained an excellent fit of the quarks mass matrices, including the CP-violating Jarlskog invariant. This leads to a unique solution for the vevs. We showed that, with the solution for the vevs obtained from the fit, it is possible to have a local minimum of the potential. We enforce this by imposing that all squared masses are positive. As in our scheme the scalar masses are not input parameters, we have to restrict one of the neutral scalars to have the mass of the known Higgs boson.

We have implemented the BFB, perturbative unitarity and the oblique parameters STU theoretical constraints. From LHC, we have considered the observed Higgs mass and the \(\kappa _V\) constraint.Footnote 5 After imposing the other constraints, we found that most of the points are close to the alignment required to respect the experimental \(\kappa _V\) constraint. We have discovered a strong correlation among the masses of the scalars, even before applying the theoretical constraints, especially for moderate to large scalar masses.

One important point is that we have numerically checked for all the points that pass our constraints, that for a given set of parameters of the potential, our minimum is the true global minimum.