Charge-breaking domain walls separating neutral vacua in multi-Higgs models

The scalar potential of a multi-Higgs model can possess a rich structure of minima and saddle points, which evolves in an intricate way as the parameters change. In the hot early Universe, it could trigger multi-step phase transitions, with exotic intermediate phases and peculiar domain wall configurations. Here, we provide a glimpse into this richness with the example of the three-Higgs-doublet model with the symmetry group $\Sigma(36)$, either exact or softly broken. We present its phase diagram tracking not only the global minimum but all of its extrema. In particular, we reveal parameter space regions in which the deepest saddle point is charge breaking. This naturally leads to phase transitions between neutral vacua which involve expanding and colliding charge-breaking bubble walls. We also comment on opportunities of multi-step phase transitions, on charge-breaking intermediate phases, and on phase transitions between different charge-breaking vacua. We illustrate the general discussion with two benchmark models, one of which possesses competing saddle points and allows for emergence of bubbles of the same true vacuum but with bubble walls of different nature. The intriguing cosmological implications of all these possibilities deserve dedicated study.


I. INTRODUCTION A. Phase transitions in multi-Higgs models
Non-minimal Higgs sectors are popular when building models beyond the Standard Model (SM), see the recent reviews [1][2][3][4].With very few assumptions, such models can address shortcomings of the SM and offer a surprisingly rich list of collider, neutrino physics, astroparticle, and cosmological signatures.A particular aspect of multi-Higgs models is a rich variety of phase transitions in early Universe they can induce, see [5,6] for recent reviews.At finite temperature T , the multi-Higgs potential can possess several competing minima, whose positions and depths evolve with T in a non-trivial way.Depending on the number of Higgs fields and the details of their interaction, this evolution can lead to multi-step phase transitions, which can be observed even within the simplistic tree-level high-T approximation of the two-Higgs-doublet model (2HDM) [7][8][9].In the past few years, multi-step phase transitions were studied in many papers within the loop-corrected effective potential, with a focus on baryogenesis and generation of non-standard spectra of gravitational waves, see for example [10][11][12][13][14][15][16][17][18] and references therein.
Phase transition sequences can include unusual intermediate phases and lead to other surprises.One example is electroweak symmetry non-restoration at high temperatures, a phenomenon known since long ago [19] and recently demonstrated within the 2HDM [17,18].Another important feature is the possibility of trapped vacua, that is, the situation when the Universe stays in a local metastable minimum all the way down to T → 0 just because the tunneling probability never gets sufficiently high to trigger bubble nucleation [20,21].
Yet another intriguing possibility is the role of charge-breaking (CB) Higgs field configurations [22,23].As pointed out in [7,8], in a certain parameter space region of the 2HDM, the Universe, while cooling down, could have evolved through an intermediate CB phase.This possibility was recently put under scrutiny in [17].Treating the finite-T one-loop effective potential of the 2HDM with the aid of the code BSMPTv2 [24], the authors found many benchmark points which exhibited charge-breaking phase at intermediate temperatures and which, at the same time, satisfied the present day collider constraints on the Higgs sector.Moreover, temperature evolution of certain benchmark points of [17] exhibited CB-to-CB first-order phase transitions within the charge-breaking phase.The very recent update of the code, BSMPTv3, confirmed the existence of intermediate CB phases in the evolution of the benchmark 2HDMs [18].
Another remarkable feature of multi-Higgs potentials is the possibility of domain walls separating spacetime domains of different minima, as well as other topologically stable configurations, which can arise from spontaneous breaking of global symmetries [25,26].If completely stable, these extended objects face severe cosmological constraints.But if they arise only temporarily, in the course of multi-step phase transitions, domain walls can play crucial role in electroweak phase transition by seeding and exponentially enhancing nucleation of the true vacuum bubbles, see for example the recent studies [27,28] and references therein.
A new twist on this topic appears when two domains of neutral minima are separated by a domain wall which is charge-breaking.The structure of domain walls linking different minima was studied through detailed numerical modeling in [29,30].It was indeed found that the path in the Higgs space providing smooth field interpolation between the two minima could contain a region with slight departure form the neutrality condition.However, the CB effect observed in [29,30] was minor.
However it is well possible to construct a 2HDM Higgs potential featuring two neutral minima separated by the deepest saddle point which is charge breaking.This is a rather straightforward exercise, and it is surprising that, for a long time, it was not discussed in literature.Only when this paper was in preparation, the publication [31] appeared, in which such CB domain walls in the 2HDM were investigated in great detail, including fermion scattering off a CB wall.
If the potential has two minima at different depths, the Universe residing initially in the false vacuum state can tunnel into the true vacuum.This transition proceeds by nucleation of bubbles of true vacuum, which then expand in the false vacuum background.The charge breaking bubble walls sweep through the primordial hot plasma and eventually collide, leading to bubble merging until the true vacuum region covers the entire Universe.Although the extensive numerical simulations of [31] offer some insights into this process, we believe that this novel regime of thermal evolution deserves a closer study, as it may offer new opportunities for baryogenesis and even generation of primordial magnetic fields [32][33][34].
Although the domain walls and intermediate phases, including charge-breaking configurations, exist within the 2HDM, we would like to check what novel features appear in even richer scalar sectors.For example, the three-Higgs-doublet model (3HDM), which was originally introduced by S. Weinberg in 1976 [35], is now being actively studied.The interest in the 3HDM is mostly driven by its richness unseen in the 2HDM, both in terms of symmetry content and model-building opportunities in the scalar and Yukawa sectors, and in terms of various phenomenological and astroparticle signals which were not compatible in the 2HDM but which coexist with three Higgs doublets, see [3] for a review.In particular, it is known that the 3HDM potential can support several minima, including coexisting neutral and CB minima [36], which may evolve with temperature in different ways.This evolution may lead, even at the tree level, to longer multi-step sequences of phase transitions with unusual intermediate phases or important roles played by domain walls, to novel opportunities for trapped vacua, and to overlapping or even competing phase transitions.All these phenomena and their possible observable consequences such as gravitational wave production and baryogenesis are worth exploring.It goes without saying that no systematic study exists of these issues within the 3HDM.

B. Goals of the present work
In this paper, we would like to stimulate this research program by considering the phase diagram and phase transitions possible within a specific version of the 3HDM based on the symmetry group Σ (36).This choice is motivated by the fact that Σ (36) is the largest finite symmetry group which can be imposed on the scalar sector of the 3HDM without leading to accidental continuous symmetries [37,38].Its scalar sector contains only four free parameters, only two of them being structurally important.This will allow us to visualize the phase diagram and to study many of its aspects analytically, including some consequences of soft breaking terms.Even in this simple setting, we can observe such phenomena as charge-breaking domain walls linking neutral minima, and CB-to-neutral and CB-to-CB phase transitions as we vary the coefficients.
We stress that the Σ(36) 3HDM is used here only as a toy model.We limit ourselves to the scalar sector and do not try to fit experimental data nor to include fermions.The main objectives of this paper are to explore the structural richness of multi-Higgs potentials beyond just looking at the global minima and to demonstrate that even in this constrained setting one find structural phenomena that have not received much attention.We believe that these features and the patterns of phase transitions are typical for multi-scalar potentials and may arise in phenomenologically viable multi-Higgs models beyond the 2HDM.We hope that the lessons learned here will guide model-builders in constructing multi-Higgs models with exotic but custom-tailored cosmological history.
The structure of this paper is as follows.In Section II we review the general properties of the Σ(36) 3HDM, describe its phase diagram, and discuss its possible global minima.Section III presents the surprisingly rich structure of extrema of the potential as a function of free parameters.We first show the results for the exact Σ(36) symmetry group.Then we add soft breaking terms which preserve some of the symmetries and explore how the structure of minima and extrema evolves as the soft breaking parameters change.In Section IV, we explore in more detail the specific issue of charge-breaking domain walls, first in the 2HDM and then in our toy model.We also give here two benchmark models which highlight two different regimes of the thermal evolution in the softly broken Σ (36).Finally, in Section V we draw our conclusions and suggest questions for follow-up studies.Additional technical results of our numerical study are presented in the Appendix.

II. Σ(36) 3HDM: GENERAL PROPERTIES
A. The scalar sector of the multi-Higgs-doublet models In the N -Higgs-doublet model, we assume that the gauge interactions and the fermion content is exactly the same as in the SM and extend only the scalar sector by adopting the idea of several Higgs generations.We introduce N Higgs doublets ϕ i , all of them having the same gauge quantum numbers, which can couple to the fermions via the standard Yukawa interactions and, importantly, interact with themselves.Their selfinteraction is described by renormalizable scalar self-interaction potential which can be generically written as all indices going from 1 to N .In the general NHDM, the parameters m 2 ij and λ ijkl are only constrained by the requirement that the potential is hermitian, which leaves room for a very large number of free parameters and renders the model hardly predictable.However, if the model is equipped with additional global symmetries forming the symmetry group G, then only a few parameters are allowed, and the model can acquire characteristic phenomenological features [3].
Let us now focus on the 3HDM.Once the Higgs potential is written, one must find its global minimum and determine the vacuum expectation value (vev) alignment of all scalar doublets ⟨ϕ i ⟩: where all vev components can, in principle, be complex.Usually, one employs the global SU (2) × U (1) transformation to set u 1 = 0 and make v 1 real, but in this paper, we will not always adhere to this procedure.
In the perturbative treatment of the model, the vev alignment defines the vacuum state; by expanding the lagrangian around the vacuum and diagonalizing the mass matrices, one obtains the physical states and their interactions.
There are two general classes of vacuum points: neutral vacuum and charge breaking vacuum.A neutral vacuum corresponds to the normal situation: the electroweak gauge group SU (2) L ×U (1) Y breaks down to the electromagnetic gauge group U (1) em , the photon remains massless, anmd the electric charge is conserved.For a neutral vacuum, there is always a global SU (2) × U (1) transformation which sets all the upper components u i = 0, so that the vev alignment becomes with v i being real or complex.If this is the case, we will often take the overall scale out, implicitly assuming that , where v = 246 GeV, and denote the vev alignment just by indicating the ratios between individual vevs.For example, one of the minima shown below (point C 1 ) has the vev alignment and label this alignment as (1, 1, 1).If the vacuum is charge-breaking, then some of the upper components u i remain non-zero in all bases.In this situation, the electroweak gauge group is broken completely, and all gauge bosons acquire masses.Although this situation does not correspond to vacuum we observe now, it could be a viable option for the early hot Universe in a finite range of temperatures, see the recent 2HDM study in [17,18].
One important remark concerns the choice of vev alignment.In the analysis of multi-Higgs models, it has become a standard practice to select a desired vev alignment as input parameters and express some of the parameters of the scalar potentials in terms if vevs.We draw the reader's attention to the fact that in multi-Higgs models with large discrete symmetry groups it is not always possible.Namely, minimization condition of such a scalar potential can only be satisfied at very special vev alignments with relations among v i , which are stable against continuous variation of the parameters of the potential.A particular consequence of this phenomenon is that a large discrete symmetry group G of the starting potential is not broken completely.In any of the possible minima, some of the original symmetries are still preserved, while the other, broken symmetries, link several distinct minima of the potential.In the case of the 3HDMs equipped with finite symmetry groups, all the patterns of symmetry breaking for all possible minima were classified in [39].

B. The potential and its parameters
The symmetry group Σ (36) was first identified as a viable option for the 3HDM in [37,38].Grouptheoretically, it is defined as a Z 4 permutation acting on the generators of the abelian group The elements of this group act on the three doublets by unitary matrices: ϕ i → g ij ϕ j , where each g can be expressed as products of powers of the generators of the group.In a suitable basis, these generators, denoted as a, b, and d, have the following expressions: where ω = exp(2πi/3).The orders of these generators are: Note that d 2 interchanges ϕ 2 and ϕ 3 ; thus, Σ(36) contains all permutations of the three doublets.Had we imposed symmetry under d 2 but not d, we would end up with the more familiar symmetry group ∆(54), first used within the 3HDMs in 1970s [40] and explored later in [41][42][43][44].It is interesting that the Σ(36) 3HDM has not yet received a detailed phenomenological study; only its main structural features are known [39,45].For more details on the relation between ∆(54) and Σ (36) and the subtleties of their definitions, see [37]; a concise version of this discussion can be found in Appendix A of Ref. [45].The 3HDM scalar potential invariant under Σ(36) has the following form: The first two lines here are invariant under all SU (3) transformations of the three Higgs doublets, and it is only the λ 3 term that selects the finite group Σ(36) out of the entire SU (3).The potential in Eq. ( 6) is also CP invariant.In fact, as was proved in [39,46], the presence of the Z 4 symmetry group within the 3HDM scalar sector automatically forbids any form of CP violation in the scalar sector, be it explicit or spontaneous.All four free parameters in Eq. ( 6) are real but play different roles.The coefficients m 2 and λ 1 only fix the overall scales of the vacuum expectation value v and the SM-like Higgs mass, but play no role in shaping the structure of the potential and in defining the number and locations of its extrema.This leaves us with only two structurally important coefficients: λ 2 and λ 3 (to be more accurate, λ 2 /λ 1 and λ 3 /λ 1 ).This is why when we study the structural features of the model we will fix m 2 and λ 1 and explore different regions of λ 2 and λ 3 .Let us also mention right away that, in order for the potential to be bounded from below (BFB), the quartic parameters must satisfy the following constraints: Note that all these inequalities are strict, which means that we impose the BFB conditions in the strong sense [69].It is known that, in a generic multi-Higgs model, one can often resort to the weaker BFB conditions on the quartic couplings which include equalities.In such a situation, the quartic part of the potential contains flat directions, but the potential is still bounded from below if the quadratic terms frow along such directions.However, in our case, the quadratic potential in Eq. ( 6) is always negative, which forces us to insist on the strong BFB conditions.

C. The global minima
An important feature of the entire ∆(54) family of the 3HDM models, including the Σ(36) 3HDM, is the "rigidity" of its global minima upon continuous variations of the parameters, a feature which is known since 1984 [47].The best way to see it is to employ the geometric method developed in [39,48].We first define , where and then rewrite the scalar potential as where In this way, we separate the "radial" and "angular" steps of the minimization problem.We first minimize Λ(x, y) ≡ Λ 1 + Λ 2 x + Λ 3 y as a function of x and y, and then, having found its minimal value Λ, we minimize the potential with respect to r 0 , which yields The SU (3) invariant quadratic term automatically leads to scalar alignment, so that we obtain a SM-like Higgs boson with the mass m 2 h = 2m 2 .The masses of the other scalars exhibit remarkable degeneracy patterns [49].
1: Global minima of the Σ(36)-symmetric 3HDM.Left: the orbit space in the (x, y) plane defined by Eqs.(8).Right: the (λ2, λ3) plane bounded by the BFB conditions in Eqs.(7), with the global minima indicated in each region (neutral and charge-breaking minima shown in blue and in pink, respectively).
What makes the form of the scalar potential in Eq. ( 9) unique is the linear dependence of V 0 on x and y.In order to find Λ, we just need to describe the orbit space, that is, the shape in the (x, y) plane formed by all possible Higgs doublet configurations, and identify its vertices.In Fig. 1, left, we show this orbit space.It has the trapezoid shape bounded by the inequalities 1/4 ≤ x ≤ 1 and 0 ≤ y ≤ x, for the derivation see Ref. [39].The value of x becomes especially clear in the bilinear formalism, which we briefly review in Appendix A and which we will use at several instances in this paper.The boundary segment corresponding to x = 1 represents the neutral vacua [50].This can be easily seen from Eq. ( 8), as a neutral vacuum point implies that all z ij = 0.The bulk of the orbit space corresponds to the charge breaking vev configurations.
Clearly, Λ must be positive everywhere in the orbit space.Writing down the condition Λ > 0 at the four vertices of the trapezoid leads to the BFB conditions listed in Eq. (7).
Minimizing a linear function of x and y in a polygon is elementary.Depending on the parameters, Λ achieves its minimum at one of the four vertices1 , labeled in Fig. 1, left.The corresponding regions on the (λ 2 , λ 3 ) plane are shown in Fig. 1, right.Below we list these four situations.
• The point (x, y) = (1, 0) corresponds to the neutral global minimum if and only if Following [39], we label a generic neutral vev alignment as (v 1 , v 2 , v 3 ) up to an overall factor: The two sets of minima B and C correspond to the same point in the orbit space of the Σ(36) 3HDM, but they represent distinct manifolds in the Higgs spaces and, therefore, different points ⃗ r is the space of bilinears, see again Appendix A. Other similar configurations such as (ω, ω2 , 1) can be reduced to the ones already listed by an overall phase rotation of the three doublets; for example, (ω, ω 2 , 1) = ω(1, ω, ω 2 ) corresponds to the alignment C 2 .Note that all six minima B and C form a single Σ(36) orbit; that is, starting form any of the six minima and applying all the transformations from the Σ(36) group, one can reach all the remaining minima.
• The point (x, y) = (1, 1) corresponds to the neutral global minimum if and only if The corresponding vev alignments are Notice that the relative phases between the vevs are not sensitive to the exact numerical values of the coefficients.It was this phase rigidity (also called calculable phases) that drove the idea of geometric CP violation in the ∆(54) 3HDM back in 1984 [47].We mention in passing that this model was analyzed in more detail in [51][52][53] and shown to be compatible with viable Yukawa sectors [54][55][56].We remind the reader that the 3HDM scalar sector invariant under the symmetry group Σ (36) does not allow one to implement any form of CP violation, neither explicit nor spontaneous [39]s.
• The point (x, y) = (1/4, 0) corresponds to the charged breaking global minimum if and only if Just as before, there are six degenerate minima in the Higgs field space.Two representative vev alignments, G and H, are Other minima can be obtained by applying the transformations from the group, just as it was done for neutral minima.
Since the three Higgs doublets transform under Σ(36) as an irreducible triplet, selecting any non-zero minimum will unavoidably break some of the symmetries.However, none of the minima breaks Σ(36) completely [39].
For each minimum point, there remains the little group of residual symmetries, which are the transformations that leave the chosen vev alignment unchanged.Further insights into the structural properties of the model, including the vev alignments, symmetry and CP properties, can be gained if one pays attention not only to the transformations from the symmetry group G but also to the transformations from SU (3) which leave G invariant, or "symmetries of symmetries" in the language of [57].The potential remains form-invariant under such transformations, only up to reparametrization of coefficients, which may provide additional links between different regimes of the same model.

III. EXTREMA OF THE Σ(36) 3HDM
A. The numerical procedure In the previous Section, we listed all vev configurations the global minimum of the Σ(36) 3HDM can have.However, we aim to track all the extrema, not only the global minima.It turns out that Mathematica cannot analytically solve the full extremization problem for the Σ(36) 3HDM potential.It becomes even more challenging for the case of the softly broken symmetry, and although useful algebraic and geometric methods have been proposed in [45,49], they have limited reach and cannot deal with the model with a generic set of soft symmetry breaking terms.
These observations force us to resort to numerical methods.For this task, we need to choose numerical values of the coefficients.Let us stress again that, in this paper, we explore structural properties of this toy model and do not perform any phenomenological fit.To this end, we set λ 1 = 1 and express all dimensional quantities in units of m 2 ; for example, the potential depth will be plotted as the reduced depths V /m 4 .Thus we only need to explore the dependence on λ 2 , λ 3 and, later, on the soft breaking parameters expressed in units of m 2 .
Once the coefficients of the potential are fixed, we find the extrema using the following numerical procedure.We parametrize the Higgs doublets via 12 real components, compute the gradient of the potential, and numerically solve the system of equations ⃗ ∇V = 0 starting from a random seed point.The solver finds an extremum, checks its neutral/charge-breaking nature, finds the eigenvalues of the Hessian, and stores this information.
We repeat this procedure N times for the same potential, each time starting from a different seed point, and compare the new extremum with the ones already stored.The value of N must be sufficiently large so that new extrema no longer appear.We found empirically that, for a typical set of the coefficients, we need up to N = 10,000 seed points to detect all the extrema.This large number of seed points looks unavoidable; limiting the search to one thousand seed points would often miss a few extrema.
More details on the our numerical procedure can be found in Appendix B.
B. 3HDM with the exact Σ(36) With the procedure described, let us now present results on the number and properties of all the extrema of the Σ(36) 3HDM potential given in Eq. ( 6).To begin with, let us fix the parameter λ 1 = 1, choose λ 2 = 1 or λ 2 = −1, and scan over the range of λ 3 from −1 to 5. The lower bound comes from the BFB conditions (7), while the upper bound is chosen to include all the λ 3 -dependent structural effects taking place in the scalar sector.
The numerical results of this scan are shown in Fig. 2 for λ 2 = 1 (left pane) and λ 2 = −1 (right pane).We show here the reduced depths of the potential V /m 4 computed for each extremum.The color indicates whether an extremum is neutral (x = 1, shown in black) or charge-breaking (1/4 ≤ x < 1, shown in shades of red).Note that when a charge breaking extremum merges with a neutral one, its color gradually shifts from bright red to black.These plots demonstrate a rich spectrum of features, which are discussed in detail in Appendix B.Here we briefly summarize the main observations and discuss how they affect phase transitions.
First, the Σ(36) potential always has many extrema, either neutral or charge breaking.As the parameters change, the depths of these points evolve, and these extrema they can merge or split.Depending on the value of λ 3 , we observed 69 or 78 extrema in total, out of which six correspond to the global minima.We remind the reader that, within the general 2HDM, the potential can admit at most two neutral minima plus four neutral and one charge-breaking saddle points, not counting the extremum at the origin.Thus, the total number of The dimensionless reduced depths of all the extrema of the Σ(36) 3HDM potentia as a function of λ3, computed for λ1 = 1, λ2 = 1 (left) and λ1 = 1, λ2 = −1 (right).The color encodes the value of x, with the neutral extrema corresponding to x = 1 (black) and the maximally charge-breaking extrema corresponding to x = 1/4 (red).
non-trivial extrema in any 2HDM is at most seven.The large number of extrema in the 3HDM comes from the higher dimensionality of the Higgs fields space.
The second observation is that the vev alignments for some extrema are rigid, in the sense that their vevs are calculable and stable upon smooth variation of λ's, while the other extrema are "floating".Among the rigid alignments, we have the global minimum, which for a positive λ 2 is attained at the points A and A ′ for −1 < λ 3 < 0 and at the points B and C for λ 3 > 0, while for a negative λ 2 , it can be either charge breaking or neutral, in accordance with the results of the geometric minimization procedure.It is also remarkable that all the special points mentioned in the previous section, even when they are not the global minima, are still present as extrema of the potential for any choice of the parameters.In addition, we found other neutral rigid points which were not identified through the geometric procedure because they never corresponded to a global minimum.A detailed discussion of the nature and regularities of these points is delegated to Appendix B.
Apart from the rigid points, we also observe in Fig. 2 other extrema, whose vev alignments and depths change smoothly as we vary λ 3 .All of these points, which we call the "floating extrema", lie in the chargebreaking part of the orbit space.Their "charge breaking depth", quantified by the parameter 1/4 ≤ x ≤ 1 and shown in Fig. 2 with the shades of red color, varies smoothly as a function of λ 3 .
We numerically confirmed that the 3HDM potential with the exact Σ(36) symmetry never develops a local minimum.We only have six degenerate global minima and saddle points.This allows us to conclude that the exact Σ(36) symmetry forbids first order phase transitions, at least at the tree level.The only phase transitions which can take place are associated with passage through points of continuous symmetries.In particular, for λ 3 = 0, the full symmetry group of the potential is promoted to SU (3), instead of the finite group Σ (36).For negative λ 2 , another transition takes place at λ 3 = λ 2 /3 = −1/3.It is also triggered by an accidental continuous symmetry because, for these values of the parameters, the phase-sensitive terms of the Higgs potential (6) can be grouped into a single combination Notice that both transitions display a large finite jump in vevs, but since the coexisting minima lie on a manifold of degenerate vacua rather than separated by a finite barrier, these transitions are not of the first order.

C. Introducing soft breaking terms
Let us now explore how the structural features of the model change if we allow for soft breaking of Σ (36).The general quadratic terms can be written as V soft = m 2 ij ϕ † i ϕ j with the hermitian matrix m 2 ij .This matrix contains nine real parameters, including the three real m 2 ii on the diagonal.A generic set of soft breaking parameters would shift the positions and depths of all the extrema, lift symmetry-induced degeneracy, and perhaps change the number of minima and stationary points.Exploring these changes in the entire 9-dimensional soft breaking parameter space and visualizing the results is a challenging task and, most likely, it will not be very revealing.
Besides, not all parameters play equal roles: some of them may trigger structural changes, while others only shift the numerical values of the observables.
Since our goal is to gain useful insights into evolution of the potential for various choices of these soft breaking parameters, we will not perform blind scans of the full 9-dimensional parameter space, but consider only such terms which preserve the Z 3 subgroup generated by the phase rotations a = diag(1, ω, ω 2 ): It is convenient to keep the overall mass scale m 2 unchanged, which implies m 2 11 + m 2 22 + m 2 33 = 0, leaving us with two independent soft breaking coefficients.We choose as free parameters the dimensionless quantities µ 1 and µ 2 , which parametrize the plane in the (m 2  11 , m 2 22 , m 2 33 ) space orthogonal to the (1, 1, 1) direction.The three coefficients m 2 ii can then be expressed in the following way: 22 33 As we will see later, this parametrization will make obvious the third-order symmetry of the (µ 1 , µ 2 ) phase diagram induced by the broken generator b in Eq. ( 5).
Let us also mention that one can organize soft breaking terms according to a different principle, proposed and developed in [49,58] and recently extended to multi-flavon models [59].Namely, one can introduce soft breaking terms which preserve the vev alignment of at least one minimum of the fully symmetric model.As shown in [49], this can be achieved if the chosen vev is an eigenvector of the matrix m 2 ij .Our choice of the soft breaking terms satisfies this condition, as all the minima of type B i are the eigenvectors of the matrix m 2 ij .However, one can in principle work with terms which softly break all the symmetries of the original model and still preserve one of the minimum; construction of such matrices was described in [49] for the same Σ(36) example as here.We delegate this more generic departure from the original symmetry to a future work.

D. The phase diagram and the evolution of extrema in softly broken Σ(36) model
Let us now explore the soft breaking phase diagram and track possible sequences of phase transitions as we move across it.Since our matrix m 2 ij is already diagonal, all three vev alignments of type B in Eq. ( 12) are its eigenvectors and, according to [49], their alignments are preserved.The values of v will, however, be different for B 1 , B 2 , and B 3 .The degeneracy among B i is lifted, and we now have three isolated local minima at different depths, separated by barriers.In contrast, the three points C in Eq. ( 13) remain at the same depth because they are linked by the unbroken symmetry a.We would like to stress here that the vev alignments C i are not the eigenvectors of the soft breaking matrix and they do evolve.By labeling an extremum as, for example, C 1 , we refer to the point which would correspond to the vev alignment (1, 1, 1) in the limit of vanishing soft breaking terms.
In Fig. 3 we present the results of a scan in the plane of the dimensionless soft breaking parameters (µ 1 , µ 2 ) for the choice of the quartic coefficients λ 1 = 1, λ 2 = 1, λ 3 = 5.The plane is composed of three "zones of influence", separated by dashed lines, each zone corresponding to the deepest minimum being B 1 , B 2 , or B 3 , as labeled on the plot.The evident three-fold symmetry of the phase diagram is due to the fact that, upon the cyclic permutation ϕ 1 → ϕ 2 → ϕ 3 → ϕ 1 , which induces the 2π/3 rotation of the diagram, the physical picture does not change, up to the permutation of the minima.
Regions with different colors correspond to models with different total numbers of minima.The central dark gray region corresponds to the potentials which still possesses all six minima B i and C i , which were present for the exact Σ(36) symmetry, although they are now located at different depths.As we move on the (µ 1 , µ 2 ) plane away from the origin, we cross a boundary, at which one of the B i minima merges with saddle points and becomes itself a saddle point.Thus, in the olive region, we have only five minima.Upon further growth of the soft breaking coefficients, three points C i cease to be the minima and, merging with low lying saddle points, become saddle points themselves.With a further increase of the soft breaking parameters, one can either stay in the light gray region, with two minima at different depths, or enter the light blue region where only one minimum survives.We checked that the boundaries at which B i successively disappear are given by straight lines, while the region where C i remain the minima is bounded by a non-trivial cycloid-like curve.
In the above diagram, we described the regions where points B i or C i remain (local) minima.When we cross their boundaries, these extrema do not disappear but just merge with other saddle points and become saddle points themselves.Since the nature of the deepest saddle point is very relevant for the phase transition dynamics, we now would like to check the evolution of all the extrema as we vary µ 1 or µ 2 .
The evolution of the reduced depths of all the extrema for the model with λ1 = 1, λ2 = 1, λ3 = 5 as a function of the dimensionless soft breaking parameters µ1 calculated at µ2 = 0 (left) and as a function of µ2 taken at µ1 = 0 (right).Black lines correspond to neutral extrema, red lines show charge breaking extrema, with the intensity of color denoting the charge breaking "depth".
In Fig. 4, we show how the depths of all the extrema change as we scan the phase diagram along µ 1 keeping at µ 2 = 0 (left plot) and along µ 2 keeping µ 1 = 0 (right plot).Starting from the Σ(36)-symmetric case and adding positive µ 1 , we observe the point B 2 to become the unique global minimum, while points B 1 , B 3 , and C to stay local minima.At the same time, the deepest saddle points, too, spread into several non-equivalent saddle points, the deepest among them approaching the minimum C.They merge around µ 1 ≈ 0.4, at the edge of the olive region in Fig. 3. Above this value, C is no longer a minimum but continues as a saddle point until about µ 1 ≈ 0.5, where it disappears as an extremum and just remains an inflection-type feature on the slope of the potential.Note that point B 1 remains a minimum until about µ 1 ≈ 0.4, even though it becomes a shallow minimum which lies higher than several charge breaking saddle points.
The right plot of Fig. 4 shows the µ 2 evolution of the extrema at µ 1 = 0; this trajectory passes vertically through the plot in Fig. 3.For this choice of soft breaking parameters, points B 1 and B 2 remain degenerate and, for positive µ 2 , they are the global minima.Note that, starting from µ 2 ≈ 0.3, the deepest saddle point is a different charge breaking point, the one shown in (18).
In general, we see that as we move from the exactly symmetric case, many previously degenerate extrema split into families of extrema at different depths.This includes, in particular, several coexisting minima, separated by barriers of different heights, so that a non-trivial sequence of first-order phase transitions can be anticipated.As the coefficients in the soft breaking coefficients grow further, many saddle points merge and then disappear, so that the structure of stationary points simplifies.This is mainly driven by the coefficients in front of some ϕ † i ϕ i becoming positive.Let us also briefly discuss the case of negative λ 3 , for which the global minimum of the symmetric model corresponds to points A, A ′ .Upon adding the soft breaking terms of Eq. ( 22), the six points A, A ′ split into a multiplet of individual minima.Although the exact vev alignments given in Eqs. ( 15), ( 16) is not preserved, we still can label them A i , A ′ i according to their fully symmetric limit.Since these soft breaking terms are CP invariant, the minima of type A and A ′ are pairwise degenerate.However each of these minima is now CP violating; thus, we observe spontaneous CP violation driven by soft breaking terms.
We also produced the phase diagram, in which, just as in the previous example, we observed several regions corresponding to different number of minima.We studied the evolution of extrema, with an equally rich list of options for phase transitions.A new feature which arises in this case is that, in large soft breaking parameters, all minima of type A and A ′ disappear and the global minimum is again attained at one of the points B. Thus, the blue regions of Fig. 3 persist even for negative λ 3 .Let us now turn to the scenario, in which the model with an exact Σ(36) symmetry exhibits a charge breaking global minimum.In Fig. 5, left, we show the phase diagram for the model with λ 1 = 1, λ 2 = −1, λ 3 = 1.Just as before, the central region, corresponding to small soft breaking terms, contains several charge breaking minima, which gradually disappear as the soft breaking terms become stronger.Away from the central region, the situation simplifies, with only one minimum remaining.It is interesting, however, that depending on the proportion between µ 1 and µ 2 , this minimum can be neutral, of type B, or charge breaking, of type G.For example, at a large positive µ 2 , the minimum can be at B 1 ∝ (1, 0, 0), B 2 ∝ (0, 1, 0), or the charge-breaking alignment

E. Charge breaking minima and phase transitions
The values ψ = 0 and ψ = π/2 correspond to the points B 1 and B 2 , respectively, while for a generic ψ, the alignment G 3 smoothly interpolates between these neutral points.This is also confirmed by Fig. 5, right, where we show the evolution of all the extrema as a function of µ 1 for fixed µ 2 = 0. We conclude that, as we shift across the phase diagram, we observe CB-to-neutral or neutral-to-CB phase transitions.They are of the second order because, as it turn out, neutral and CB minima do not coexist in this example.They can become the first order beyond the tree level, as it happened in the 2HDM when the one-loop effective potential treatment of [17] improved over the tree-level approach of [7,8].However we do observe the first order phase transitions inside the CB phase even at the tree level, just because several CB minima can certainly coexist at small soft breaking parameters.
A general take away message is that all these structural opportunities, with their possibly intriguing phase transitions dynamics, arise thanks to soft breaking terms.

IV. CHARGE BREAKING BUBBLE WALL BETWEEN TWO NEUTRAL MINIMA A. Charge breaking bubble wall in the 2HDM
As we observed in Section III, the Σ(36) 3HDM potential with λ 2 > 0 and a sufficiently large λ 3 contains the deepest saddle point which is charge breaking, not neutral.Before exploring it, we would like to remark that such a situation does not require three Higgs doublets and is possible within the 2HDM.This simple fact remained for a long time unexplored and, to our best knowledge, it was first explicitly pointed out only in the very recent publication [31].To give a concrete example tractable analytically, let us consider the 2HDM potential It is invariant under the sign flips of each doublet, ϕ 1 ↔ ϕ 2 , and the usual CP transformation, which combine, group theoretically, into the (Z 2 ) 3 group acting in the space of bilinears [60,61].This potential can also be obtained by imposing a single generalized CP symmetry [62][63][64].Let us also assume that the quartic parameters λ i satisfy, in addition to the usual BFB conditions, the following inequalities: λ 3 > λ 1 and λ 4 − |λ 5 | > 0.Then, straightforward algebra allows us to find all the extrema of this potential.We get a pair of degenerate global minima at two pairs of neutral saddle points at as well as a charge breaking saddle point at With the parameters satisfying λ 3 > λ 1 and λ 4 − |λ 5 | > 0, the charge breaking extremum becomes the deepest saddle point.By adding soft breaking terms, we can avoid degeneracy between the two minima while keeping the CB saddle point the deepest one.If the early Universe temporarily resided in the local metastable minimum, it could decay, via tunneling or thermal fluctuations, to the true minimum, see for example the pedagogical expositions in [5,65].As a result, bubbles of true vacuum will spontaneously form and expand in the false vacuum background, sweeping across the primordial hot plasma and eventually colliding and merging.Not aiming at any accurate computation of the nucleation rate, we only mention that, neglecting thermal fluctuations and treating the tunneling integral in the thin-wall approximation, we find the exponentially suppressed value exp(−S E ), where S E is the euclidean action calculated for the so-called bounce solution [66].φ V δ ϵ FIG.6: A schematic profile of a titled double-well potential.
In the simplest example of one real scalar field with a tilted double-well self-interaction potential with the quartic coupling λ as schematically depicted in Fig. 6, S E depends on the height of the barrier δ and the difference ϵ between the depths of the two minima, in the following way: Here, C is the numerical coefficient of the order of a few hundred, which is inessential for the present qualitative discussion.A good rule of thumb is that, for reasonable values of λ, the values δ/ϵ ≪ 1 lead to a very quick tunneling, while δ/ϵ ≫ 1 are characteristic of extremely long-lived trapped metastable states.With these preliminary remarks, let us now consider a multi-field potential, with two different minima and several competing saddle points.It may happen that several paths passing through different saddle points exist.For each competing path, the action involves the barrier height δ of its saddle point and an effective λ, which essentially depends on how long this optimal path is.Although these paths are not necessarily straight, we can roughly estimate 1/λ ∼ (∆v) 4 /δ, where ∆v is the straight segment distance between the minimum and the saddle point.These heuristic arguments lead us to the general conclusion that, among several competing paths, the minimal action is achieved along the path passing through the deepest saddle point, provided it does not lie too far from the two minima it connects.Applying this qualitative criterion to the 2HDM example given above, we see that the path passing through the charge breaking saddle point is indeed optimal.
Thus, the phase transition will proceed through formation, expansion, and coalescence of bubbles of true vacuum with CB bubble walls.It is interesting to see whether these CB bubble walls, sweeping through plasma, has any consequences for baryogenesis or leads to fermion charge separation, capable of seeding large-scale charge inhomogeneities and leading to primordial magnetic fields.The recent study [31] provided only a few first insights into possible phenomena which could take place during this disruptive process, and more studies are definitely needed.

B. Charge breaking bubble walls in the Σ(36) 3HDM
Returning to the Σ(36) 3HDM, let us first remark that we have identified two regimes where neutral minima of different depths are separated by a charge breaking domain wall.The first one exists already for the exactly symmetric case, while the second one appears upon adding a significant soft breaking term.
First, when studying the phase diagram for the softly broken case, we selected λ 3 = 5 which guaranteed that the deepest saddle point was charge breaking, see Fig. 2.This is a floating point, so that its vev alignment smoothly depends on the values of the coefficients, but in the limit of large λ 3 , it asymptotically approaches This asymptotic vev alignment corresponds to the point (x, y) = (1/2, 0).Of course, by applying Σ(36) transformations, we can find other vev configurations corresponding to this point; in total, there are nine.
Let us now add soft breaking terms with moderately large µ 1 passing from negative to positive values, while keeping µ 2 = 0. Initially, the global minimum is at B 1 .As we pass the symmetric point and move into the positive µ 1 region, B 1 becomes a local minimum ready to tunnel to a lower-lying minimum, either C, B 3 , or the global minimum B 2 .To which minimum it will actually tunnel is not immediately clear and requires a dedicated numerical study.What is certain is that the saddle point through which this transition occurs is charge breaking, and the natural expectation is that it is the point in Eq. ( 29).
If µ 1 increases further, this deepest saddle point approaches the neutral minima of type C, merges with them resulting in neutral saddle points, see again Fig. 4, left.However if, in addition to the moderately large µ 1 , we add a sufficiently large positive µ 2 , we enter the second regime with CB saddle points, similar to Fig. 4, right.This time, the two competing minima are of type B 1 and B 2 , while the CB saddle point is of type G 3 , as in Eq. ( 23).This regime, which closely resembles the 2HDM example shown above, persists even for large µ 2 , and the deepest CB saddle point never merges with the neutral vev manifold.

C. Benchmark models
To illustrate the above general discussion, we provide in this subsection two benchmark models with softly broken Σ (36), which lead to first order phase transitions involving charge breaking bubble walls.In the first model the CB nature of the wall in unambiguous, while in the second model we observe several competing saddle points of different nature, which potentially leads to the remarkable situation of several bubble of the same new vacuum separated by different bubble walls.
We begin by writing, at zero temperature, the quadratic part of the potential as where the coefficients m2 ii ≡ −m 2 + m 2 ii include the overall m 2 from the fully symmetric potential Eq. ( 6) and the soft breaking terms of Eq. (21).We adjust the values of m 2 and λ 1 to reproduce v = 246 GeV and m h = 125 GeV.As suggested by the analysis in the previous section, we often encounter minima of type B i .Their vevs are the eigenvectors of the soft breaking terms we use, which leads to scalar alignment even in the softly broken case [45].This allows us to fix the parameters m 2 = m 2 h /2 and λ 1 = m 2 /v 2 ≈ 0.129 at zero temperature.Just as in Section III D, we also fix λ 2 = λ 1 and λ 3 = 5λ 1 , so that we can rely on the phase diagram shown in Fig. 3 and the extrema evolution in Fig. 4.
The two benchmark models then differ only by the zero-temperature values of the soft-breaking parameters µ 1 and µ 2 : The two choices of (µ 1 , µ 2 ) are indicated in Fig. 3, by circles with arrows.The first model corresponds to the point just below the right dashed line, while the second model refers to the point near the boundary between the black and blue regions.In both cases, the zero-temperature minimum is at the vev alignment B 3 = (0, 0, 1).We now add temperature evolution.At non-zero temperature, the shape of the potential changes, which can be analyzed within the one-loop thermal effective potential [67,68].To make the physical consequences transparent, let us focus on the main effect in the high-T approximation, the thermal T 2 contributions to the scalar masses.These are induced by one loop corrections to scalar field propagators: Once we know m2 ii and c i , we can represent the thermal history of the scalar sector as a trajectory on the (µ 1 , µ 2 ) plane, see a 2HDM example in the recent study [17].These trajectories are straight, and their directions for the two benchmark models are shown by the arrows in the Fig. 3.
Due to the Σ(36) symmetry of the quartic interactions, the scalar one-loop contributions are identical for all three Higgs doublets.The gauge boson contributions are also identical and equal to the standard expression for electroweak doublets.Thus the coefficients c i can only differ in how each doublet couples to fermions: We explicitly indicated here the color factor N c = 3 for quarks and N c = 1 for leptons.The fermion loop contributions depend on the Yukawa structures and, in particular, on the top-quark coupling with the three Higgs doublets.In both benchmark models, we assume that the top-quark couples only to the third Higgs doublet with the SM-like coupling: y Of course, other Yukawa sector models are possible and should be included when building a realistic Σ(36) 3HDM model.In Fig. 7, we show temperature evolution of the depths of the minima B 2 and B 3 as well as the few lowest lying saddle points in the benchmark model 1 (left) and model 2 (right).For model 1, we show only the deepest saddle point, as all additional saddle points lie significantly higher.Let us track the evolution of the Higgs potential as the Universe cools down.At very high temperatures, the quadratic coefficients are all positive and the electroweak symmetry is restored, the minimum staying at v = 0. Around T = 120 GeV, the global minimum shifts to the alignment B 2 = v(T )(0, 1, 0), with the value of v(T ) gradually increasing as the temperature drops.Below T ≈ 90 GeV, the vev alignment B 3 transforms from a saddle point to a minimum, and its depth quickly grows.At 70 GeV, B 3 becomes the global minimum.However the Universe still resides in the local minimum B 2 due to a significant potential barrier which separates the two minima.This potential barrier goes over the deepest saddle point which is charge breaking and is shown in Fig. 7 in red.
The actual phase transition happens at a lower temperature, when the difference of the depths of the two minima becomes comparable to the saddle point barrier height, or δ/ϵ ∼ 1 in the notation of Section IV A. To guide the eye, the pink dashed line in Fig. 7, left, indicates the temperature when δ/ϵ = 1, which is around 55 GeV.Note that we do not claim that the phase transition occurs exactly at this point, for such a claim would require an exact numerical computation of the bounce integral in the multi-Higgs space and the evaluation of the bubble nucleation rate.We just point that it will happen at temperatures somewhat below the point of the B 2 /B 3 crossing, and that it will be a first order phase transition.Since the deepest saddle point is charge breaking and all other saddle points lie significantly higher, this phase transition definitely proceeds via emergence of bubbles of the vacuum B 3 in the background of the false vacuum B 2 , separated by charge breaking bubble walls.This is exactly the regime predicted by our general analysis.
In the second benchmark model, shown in Fig. 7, right, temperature evolution proceeds similarly.However, the choice of the zero-temperature point (µ 1 , µ 2 ) now leads to a much more intricate evolution of the saddle points.We see that in the temperature region where a phase transition is expected to occur we have two families of saddle points, plus yet another family of local minima of type C. We find that the two competing saddle points are charge breaking, although one of them is close to a neutral Higgs configuration, and are well separated in the Higgs field space.These extrema display similar but not exactly identical temperature evolution.Thus, within the benchmark model 2, our simplified analysis is unable to answer which of the two competing saddle points gives the higher nucleation probability.It is, in fact, well possible that new vacuum regions with both types of bubble walls emerge, expand and merge.
The lesson we learn from benchmark model 2 is that situations with competing saddle points are rather generic in the symmetry-based 3HDMs.Our simplistic analysis only indicates that such peculiar regime indeed exists but cannot accurately predict the dynamics of phase transition.Quantitative understanding of the evolution of the Universe, as well as its implications for gravitational wave productions, baryogenesis, and primordial magnetic field generation, requires a dedicated numerical study at the level of finite-temperature one-loop effective potential.

V. DISCUSSION AND CONCLUSIONS
In this work, we presented an exploratory study of a multi-Higgs potential with a large symmetry group, either exact or softly broken, in which we not only looked at the global minimum but also tracked all extrema of the potential.The calculations were done for Σ(36) 3HDM, but we believe the lessons learned can be useful in other multi-Higgs models.Here are the main observations made in the course of this study.
• Potentials with large finite symmetry groups, even with very few free parameters, possess surprisingly rich structure of minima and saddle points.In the fully symmetric 3HDM case, we observed up to 78 extrema, either neutral and charge breaking, which is a significant increase with respect to the 2HDM, where the potential can possess at most seven stationary points.Vev alignments of the extrema can be "rigid" (stable against smooth variation of the free parameters) or "floating" (smoothly shifting as the parameters change).There are also bifurcation points where saddle points merge or split, with interesting symmetry driven patterns.
• When we searched for all the stationary points via numerical methods, we were forced to use several thousand seed points in order not to miss any extremum.This was a significant bottleneck in our study, and it calls upon more efficient methods for locating all stationary points of a multi-Higgs potential.
• There are regions in the parameter space where the deepest saddle point is charge breaking.As a result, there can exist charge-breaking domain walls separating neutral minima.We also gave an explicit 2HDM example with the same feature.
• In the softly broken case, we found two distinct regimes, in which tunneling between two neutral minima at different depths proceeds via expansion and collision of charge-breaking bubble walls.
In short, we find it remarkable that such a structurally simple model as Σ(36) 3HDM with soft symmetry breaking allows for several phase transition regimes without any fine-tuning.We believe that a detailed numerical investigation of many of these processes may lead to new insights into possible dynamics of the early Universe.
A particularly intriguing phenomenon is tunneling between neutral vacua which proceeds via expansion and collision of charge-breaking bubble walls.Although this feature was recently identified in [31] within the 2HDM, the multi-Higgs example we considered here offers novel opportunities.First, we observe two different regimes for this phenomenon, corresponding to weak and to strong soft breaking terms, with potentially different consequences.Second, the 3HDM example possesses many more minima and saddle points than the 2HDM, with several competing phase transitions and even competing tunneling paths between the same pair of minima.
To illustrate some of these findings, we presented two benchmark models, including one in which saddle points of different nature indeed compete.This situation can potentially lead to a remarkable cosmological scenario in which bubble of the same true vacuum appear and expand in the false vacuum background, but these bubbles possess bubble walls of different nature.
Numerous questions for follow-up studies naturally emerge.Staying within the same toy model, the softly broken Σ(36) 3HDM, one can explore in more detail the exact sequence of phase transitions.In particular, for a potential with generic soft breaking terms, one needs to identify the exact trajectory in the Higgs field space which minimizes the euclidean action between a chosen pair of minima.As several saddle points can compete, it may happen that the bubble nucleation rate could depend in a highly non-trivial way on the parameters of the potential.Also, if the symmetry group is not fully broken, one can perhaps observe the situation when different bubbles produce different true vacua at the same depth.It will be interesting to see whether genuine stable domain walls can be produced or whether one of the competing vacua eventually takes over.One can also imagine the situation with two phase transitions taking place in rapid succession, so that the second one sets in before the bubbles from the first phase transition have already expanded and fully merged.Extensive numerical simulations are required to verify all these peculiar evolution regimes.
Next, one can aim at building a realistic model incorporating the above phase transition dynamics.Using the results of this paper as guidelines, one can construct a 2HDM or a 3HDM with a softly broken symmetry group, complete it with a suitable Yukawa sector, and use the finite-temperature loop-corrected effective potential to check whether the regimes outlined here persist in the actual thermal history of the model.Multi-step phase transitions, competing minima and saddle points, dynamics of charge-breaking bubble walls, charge-breaking, charge-restoring, or CB-to-CB phase transitions are the features to look for.The recent study [17] of the charge breaking vacuum at intermediate temperatures with all sorts of accompanying phase transitions proves than the one-loop effective potential picture can be even richer than the tree-level one.
Finally, one should explore in depth whether the tantalizing possibility of charge-breaking bubble walls and charge-breaking intermediate phases has significant impact on baryogenesis, on gravitational wave spectrum, and on the origin of primordial magnetic fields [32][33][34].It may well happen that the intricate multi-Higgs dynamics discussed here opens novel opportunities for the cosmological history of the early Universe.
• As we change the coefficients of the potential, or add soft breaking terms, the extrema can shift, merge, or split.The depth of the potential at each point also changes.We found it necessary to track the evolution of all the extrema in order to have the good understanding of the global picture.
b. Rigid vs. floating extrema In addition to the main text, we provide here more details on the non-trivial behavior of the extrema in the Σ(36)-symmetric 3HDM.To this end, we reproduce in Fig. 8 the left plot of Fig. 2 together with a zoom on the narrow range of V /m 4 where mergers and splittings of various extrema take place.
The first observation is that the vev alignments for some extrema are rigid (their vevs are calculable and are stable upon smooth variation of λ 3 ), while the others are "floating".Among the rigid alignments, we have the global minimum, which is attained at points A and A ′ for −1 < λ 3 < 0 and at points B and C for λ 3 > 0. It turns out that all the special points mentioned in the previous section, both neutral (A, A ′ , B, C) and charge-breaking ones (F, F ′ , G, H), even when they are not minima, always represent rigid extrema of this potential.Their vev alignments are not sensitive to the exact values of the coefficients, and the depths of the potential V /m 4 at these points follow simple laws: In addition, we found other neutral rigid points which were not identified through the geometric procedure of the previous section because they lie inside the orbit space and never correspond to a global minimum.The Other points of the same type can be obtained by applying Σ(36) transformations.Notice that the extremum s − has the value y ≈ 0.933, which is remarkably close to 1.This explains why this curve follows so closely the vev alignment AA ′ , which corresponds to y = 1.
In addition to the rigid points, we also observe in Fig. 2 other extrema, whose alignments and positions change smoothly as we vary λ 3 .All of these points, which we call the "floating extrema", lie in the chargebreaking part of the orbit space.Their "charge breaking depth", quantified by the parameter 1/4 ≤ x ≤ 1 and shown in Fig. 2 with the shades of red color, varies smoothly as a function of λ 3 .
Since we keep track, for each extremum, of the eigenvalues of the Hessian, we observe that, for positive λ 2 , only the previously found points A, A ′ , B, C can be the minima.There are no other local minima in this potential; all other extrema are saddle points, although with different signatures of the Hessian matrices.This allows us to conclude that within the 3HDM with the exact Σ(36) symmetry, no first order phase transitions are possible, at least at the tree level.
The only phase transition which can take place for positive λ 2 is the transition from A, A ′ to B, C at the point λ 3 = 0, the value at which the Σ(36) symmetry is accidentally restored to the full SU (3).Although there is a large finite jump in vevs as we move through this point, this is not a first order phase transition since the coexisting minima are not separated by a finite barrier.

c. Bifurcation points and the number of extrema
Unlike rigid extrema, which exist for all values of the parameters, the floating extrema are present only within certain ranges of λ i .As we scan over λ 3 , we observe them either merge with or emerge from a rigid neutral extremum.A closer inspection shows that these bifurcations happen only along the curves corresponding to the rigid saddle points s 0 , s + , s − , not the global minima.We do not have the analytic expressions for the floating point location, but we do have them for the rigid points.In order to find the exact instants of splitting or merging, we use the condition that, at the bifurcation point of the curves of s 0 , s + , s − , one of the charged Higgs masses vanishes.The charged Higgs mass matrix can be computed analytically, thus we set its determinant to zero and obtain an equation for λ 3 .Solutions of these equations give the values of λ 3 at which

FIG. 3 :
FIG.3:The phase diagram of the softly broken Σ(36) model in the plane of dimensionless soft breaking parameters (µ1, µ2) for λ1 = 1, λ2 = 1, λ3 = 5.The color encodes the total number of minima; when the minimum is unique, it is explicitly labeled.The two circles with arrows correspond to the thermal evolution of the two benchmark models discussed in Section IV C: benchmark model 1 (blue) and benchmark model 2 (red). µ

2 µ 2 = −0. 2 B 2 B 3 FIG. 7 :
FIG.7: Temperature dependence of the depths of the minima (shown in blue) and the few deepest saddle points (neutral shown in black, charge-breaing shown in red) of the Higgs potential in the two benchmark models, see the main text for the parameters.The vertical pink dashed lines show the temperatures when the barrier height due to the deepest saddle point is equal to the difference between the depths of the two minima.

FIG. 8 :
FIG.8:The reduced depths of all the extrema of the Σ(36) 3HDM as a function of λ3, computed for λ1 = 1, λ2 = 1.Left: the full picture.Right: a zoom on a narrow region of V /m 4 .Color encodes the value of x, with the neutral extrema corresponding to x = 1 (black) and the maximally charge-breaking extrema corresponding to x = 1/4 (red).The labels indicate some of the rigid extrema, see the main text.