Domain Wall Constraints on Two Higgs Doublet Models with $Z_2$ Symmetry

The Two Higgs Doublet Model (2HDM) with spontaneously broken $Z_2$ symmetry predicts a production of domain walls at the electroweak scale. We derive cosmological constraints on model parameters for both Type-I and Type-II 2HDMs from the requirement that domain walls do not dominate the Universe by the present day. For Type-I 2HDMs, we deduce the lower bound on the key parameter $\tan\beta>10^5$ for a wide range of Higgs-boson masses $\sim$ 100 GeV or greater close to the Standard Model alignment limit. In addition, we perform numerical simulations of the 2HDM with an approximate as well as an exact $Z_2$ symmetry but biased initial conditions. In both cases, we find that domain wall networks are unstable and, hence, do not survive at late times. The domain walls experience an exponential suppression of scaling in these models which can help ameliorate the stringent constraints found in the case of an exact discrete symmetry. For a 2HDM with softly-broken $Z_2$ symmetry, we relate the size of this exponential suppression to the soft-breaking bilinear parameter $m_{12}$ allowing limits to be placed on this parameter of order $\mu$eV, such that domain wall domination can be avoided. In particular, for Type-II 2HDMs, we obtain a corresponding lower limit on the CP-odd phase $\theta$ generated by QCD instantons, $\theta \ \stackrel{>}{{}_\sim}\ 10^{-11}/(\sin\beta \cos\beta)$, which is in some tension with the upper limit of $\theta \ \stackrel{<}{{}_\sim}\ 10^{-11}$--$10^{-10}$, as derived from the non-observation of a non-zero neutron electric dipole moment. For a $Z_2$-symmetric 2HDM with biased initial conditions, we are able to relate the size of the exponential suppression to a biasing parameter $\varepsilon$ so as to avoid domain wall domination.


Introduction
Domain walls are topological defects which emerge from the breaking of discrete symmetries [1] resulting in a vacuum manifold containing topologically disconnected points. These disconnected points correspond to degenerate vacua. Phase transitions producing domain walls occur at a finite rate and as such, the field can select different vacua in causally disconnected regions of space. This divides the universe into so-called domains the interfaces between which are domain walls [2]. Defects which emerge from the breaking of a global symmetry are expected to enter a regime of dynamical scaling such that the number of defects is constant per Hubble horizon [3]. Domain walls follow a power law scaling with an exponent close to 1 as shown in simulations for the so-called Goldstone model with a single real scalar field [1,4]. This is due to the fact that the walls have a tension under which they collapse as quickly as causality permits. This scaling feature results in an undesirable fate for the Universe where domain walls can be present in nature, in stark contrast to all our observations. The energy density of matter and radiation both scale proportionally to (time) −2 in their respective epochs of domination. However, domain wall energy density scales proportionally to (time) −1 [5]. This means that domain walls will come to dominate the Universe at late times [6][7][8]. This is the so-called domain wall problem. Therefore, if cosmic domain walls are to exist in nature constraints must be placed on domain wallforming models such that domain wall domination does not occur [5] or, at least, occurs after present day. Domination could be avoided by, for example, having an additional field couple to the walls altering their scaling behaviour [4,9,10]. Alternatively, one could have the domain walls decay before they come to dominate the Universe by making the discrete symmetry approximate. It is important to note that these modifications require a change in the symmetries of the model and must, therefore, be well-motivated given the fundamental role of symmetries in physics.
The Higgs mechanism for electroweak symmetry breaking was verified by the measurement of a Higgs boson [11,12] of mass 125 GeV at the LHC [13]. The properties of this scalar particle so far match those predicted for the SM Higgs scalar [14,15]. Nonetheless, current experimental measurements do not prohibit the existence of more scalar particles. One minimal and theoretically well-motivated extension which can be made to the SM is to introduce a second complex Higgs doublet into the theory. This is the so-called Two Higgs Doublet Model (2HDM) first suggested by T. D. Lee in 1973 [16]. A complete study of the general CP-violating 2HDM is given in [17] (for a review, see [18]).
It is well known that the 2HDM predicts the emergence of a variety of topological defects, such as domain walls, vortices and global monopoles, from the breaking of accidental symmetries which the model can posses under certain parameter choices [19][20][21][22]. In the thermal history of the Universe, theories of New Physics, such as the 2HDM, can predict a series of symmetry breaking phase transitions as the Universe expanded and cooled [2]. These broken symmetries are no longer observable but should be restored in the early Universe when temperatures were far higher than at present [1]. These phase transitions can leave relic topological defects which can serve as probes of high energy physics in the early Universe [1,3,23].
In this article we focus our attention on the 2HDM with Z 2 symmetry, whose spontaneous breaking predicts the existence of domain wall solutions. In particular, we obtain cosmological constraints on the theoretical parameters of the 2HDM which arise from the non-observation of such domain walls. By solving the relevant equations of motion [1,4,9,24], we present a number of numerical simulations of such topological defects that consolidate the cosmological constraints derived in this paper.
The remainder of this article is structured as follows. In Section 2, we outline the scalar sector of the 2HDM with softly-broken Z 2 symmetry, introduce the physical degrees of freedom in the model and provide a brief review of experimental constraints on the 2HDM for models of Type I and II. In Section 3, we present constraints on the 2HDM from domain wall domination for cases where the Z 2 symmetry is exact and parameter regimes in which domain wall domination could be avoided. In Sections 4 and 5, we present results of simulations of the 2HDM with both an approximate Z 2 symmetry and biased initial conditions for a 2HDM where the Z 2 symmetry is exact. In both cases, domain wall networks are unstable and domain wall domination can be avoided by requiring that these domain wall networks be sufficiently short-lived. Finally, Section 6 summarizes the main results of our study.
The vacua are parametrized as This parametrization will be referred to as CP-preserving vacua. The parameters v 1 , v 2 , v + and ξ are referred to as the vacuum manifold parameters. For the CP-preserving vacua (2.3), the VEVs can be calculated in terms of the potential parameters, where we have definedλ 345 = λ 3 + λ 4 − |λ 5 |. The 2HDM has 5 physical scalar particles: 2 neutral CP-even states, h and H, one CP-odd neutral state, A, and 2 charged states, H ± . The other three scalar degrees of freedom correspond to would-be Goldstone bosons, G 0 and G ± , which are absorbed into the longitudinal components of of the electroweak gauge bosons, W ± and Z 0 . We identify h with the Higgs particle measured at the LHC by ATLAS and CMS [11,12] fixing the parameter M h at 125 GeV [13]. Furthermore, the SM VEV is fixed at v SM = 246 GeV.
In order to investigate the evolution of domain walls in the 2HDM with approximate Z 2 symmetry we first obtain a physical parametrization of the model with which to perform our numerical simulations. Expressions for the masses of the scalar Higgs particles, h, H, A and H ± , are obtain as eigenvalues of the Hessian matrix of (2.2) using the parametrization where ϕ + i are complex scalar fields. The CP-even mass matrix is derived to be while the CP-odd mass is where we have introduced the short-hand notations sin x = s x and cos x = c x . The charged Higgs mass can then be written as The CP-even mass matrix is diagonalized by the mixing angle, α, (2.9) Using the above expressions one obtains the physical parametrization of the scalar potential Rescaling for dimensionless energy per unit area,Ê = E/M h v 2 SM , we can write (2.2) in a dimensionless form, , and α/β are the CP-even/odd mixing angles which diagonalize the scalar mass matrices. For details of this reparametrization and rescaling procedure, see [10].
In order to have a phenomenologically acceptable model, we must also consider that current measurements of signal rates for the Higgs discovered at the LHC are close to those predicted by the SM [14,15]. Therefore, we must restrict our investigation to parameters in/near the so-called SM alignment limit where the couplings of the CP-even scalar, h, match those predicted by the SM. Exact SM alignment is obtained when the relation cos (α − β) = 1 holds 1 . Therefore, the physical parametrization we will use in phenomenological discussions to follow will be . When choosing physical parameters and considering the implications of our phenomenology in later sections we must also account for the effect of model type on current experimental constraints on the 2HDM. The Type of a 2HDM is determined by the form of its Yukawa sector. In other words, constraints on the masses of the 2HDM scalars from experimentally-measured signal rates are type-dependent due to differences in the Higgsfermion interactions of the models. The Yukawa Lagrangian can be written in its most general form as where q L and L are SU(2) L doublets for left-handed quarks and leptons, repsectively; u R , d R and R are SU(2) L singlets for right-handed up-type quarks, down-type quarks and leptons, respectively 2 and Y u,d, i are 3 × 3 Yukawa matrices for each Higgs doublet. However, the 2HDM Yukawa sector as given in (2.12) is too general and restrictions on the Higgs-fermion couplings are required to remove or limit tree-level flavour-changing neutral currents (FCNC) [18]. The two models we will consider are the so-called Type I and Type II 2HDMs [29]. In Type I models all fermions couple to only one of the doublets (conventionally chosen to be Φ 2 ) [18] whereas in Type II models, up-type quarks couple to one doublet (Φ 2 by convention) and down-type quarks and leptons to the other [25] 3 .
For a 2HDM with softly broken Z 2 symmetry of Type I 2HDM, flavour physics constraints place a limit of tan β > 1 for M H± = 1 TeV with the constraint strengthening to tan β > 3 for M H ± = 100 GeV [25,26]. As such, the entire range of masses from 100 GeV upwards can be chosen from M H ± without contradicting flavour constraints provided sufficiently large values for tan β are chosen. Constraints on the charged Higgs from direct and indirect detection at the LHC increase the lower bound on tan β for M H ± < 300 GeV [25]. For a Type II 2HDM flavour physics constraints place much stronger bounds on the charged Higgs mass. Specifically, a lower bound of M H ± 600 GeV exists for all tan β 1 Constraints from experimental limits on the mixing angles for varying values of the Higgs masses can be found in [25][26][27][28]. 2 Note that all these objects are 3-vectors is flavour space where flavour indices have been suppressed. 3 The specific forms of the Higgs-fermion couplings in each case can be found in Table. 2 of [18] in terms of the mixing angles, α and β.
and M H ± 650 GeV for tan β < 1 independent of the other physical 2HDM parameters [25]. In both Type I and Type II 2HDMs the combined constraints of [26] require a strong alignment between M H ± and either M H or M A . This is attributed to the strong constraint on the value of cos(α − β) from current Higgs coupling measurements being close to SM alignment. As with the charged Higgs constraints, the Type II model is more strongly constrained by current observations with the entire parameter range of the neutral Higgs masses, 100 GeV ≤ M H,A ≤ 1000 GeV considered in [26], ruled out for their lower benchmark charged Higgs masses of M H ± = 250 GeV and 500 GeV. Therefore, we choose alignment of the charged Higgs mass with either the scalar mass, M H , or the pseudoscalar mass, M A , motivated by the combined constraints of [26].

Exact Z 2 Symmetry
In an FRW universe, the energy density of both matter and radiation decrease proportionally to (time) −2 in their respective epochs of domination. However, domain wall energy density decreases proportionally to (time) −1 . Therefore, the energy density of domain walls will increase relative to matter and radiation in their respective epochs and, hence, come to dominate the energy density of the universe. The time of this domination is determined by the energy per unit area of the domain walls. It is clear that we do not live in a domain wall dominated universe and, therefore, any model which produces domain walls must not allow domination before present day. It has been established that domain wall networks in the Z 2 -symmetric 2HDM exhibit a deviation from ∝ t −1 scaling [10]. Specifically, more domain walls are predicted at late times in the 2HDM than one would expect for ∝ t −1 scaling. This feature makes the domain wall problem more restrictive. Here, we calculate the domain wall density for the Z 2 -symmetric 2HDM assuming ∝ t −1 scaling and require that domain wall domination occurs after present day to obtain corresponding constraints on the physical observables. As such, this calculation provides a minimal constraint from 2HDM domain wall domination.
In the Z 2 -symmetric 2HDM, the energy per unit area of the domain walls is given by for CP-preserving vacua (2.3). The topologically non-trivial solution which minimizes the energy per unit area can be obtained via gradient flow (see, for example [10,19]). It should be noted that the SM VEV, v SM = v 2 1 + v 2 2 , also changes in the vicinity of the kink as the solution interpolates from one vacuum to another. The energy density of a domain wall network can be approximated by a self-scaling argument. The total energy within a Hubble horizon of radius, r is proportional to Er 2 . Therefore, the energy density, ρ dw ∝ Er −1 and since the horizon expands at the speed of light, ρ dw ∝ Et −1 . Hence, we write where A is a constant of proportionality, quantifying the number of walls per horizon, and define a corresponding domain wall density parameter in the usual manner, with critical density at present day where we have used H 0 = 72km s −1 Mpc −1 = 1.54 × 10 −42 GeV in natural units. For Ω dw < 1 at present day, we obtain the limit Therefore, for t 0 6.6×10 41 GeV −1 and M pl 1.2×10 19 GeV, we obtain the dimensionless inequality, Firstly, it should be noted that agreement with this limit can always be obtained for sufficiently large or small values of tan β. It is always energetically favourable for the kink to interpolate between the smaller of the two VEVs. Since tan β determines the relative size of the VEVs, for tan β > 1 the first doublet is Z 2 odd while for tan β < 1 the second doublet is Z 2 odd. This is illustrated in Fig. 1 for some benchmark values of the CP-even scalar mass in the SM alignment limit. We see in the left panel of Fig. 1 that domain walls do indeed become ultra-light in large and small limits of tan β where the VEV of the Z 2 odd doublet becomes vanishingly small. The parameter tan β is the primary parameter in the variation ofÊ. In the right-hand panel of Fig. 1 we see that the alignment parameter cos(α − β) only has a weak effect on the energy density. Moreover, the impact of changing M H is also weak.
The variation of dimensionless energy per unit area,Ê in SM alignment is given in Fig. 2. Note that there is some subtlety in the impact of the parameter A on this limit. This is more than simply an assumed number of domain walls per horizon as this value will change between the matter and radiation eras through which these domain walls scale. Nonetheless, this should note affect the order of magnitude of the limit (3.7). We find that for domain walls of the type seen in the neutral vacuum minimum energy kink solutions a domain wall problem can only be avoided for experimentally viable Higgs masses at large values of tan β (of the order 10 5 or more). In lower tan β regimes one cannot evade the constraints placed on the Z 2 -symmetric 2HDM by domain wall domination without unreasonably low values of the scalar masses. It should be made clear that these results only pertain to scenarios where the 2HDM possess an exact discrete symmetry and, hence, a domain wall problem. These stringent constraints suggest that in order to have cosmologically viable 2HDM domain walls in experimentally viable parameter regimes a means of modifying the scaling behaviour of these domain walls will be required.

Approximate Z 2 Symmetry
So far, we have only considered the scenario in which the 2HDM possesses an exact Z 2 symmetry. We have shown that the domain wall problem present in this model places highly restrictive limits on the parameters of the model such that domination can be made to occur after present day. We now turn our attention to means of eliminating the domain wall problem altogether. For the 2HDM with softly-broken Z 2 symmetry, (2.2), the degeneracy of the vacua is removed and the scalar potential contains so-called true and false vacua. The true vacuum is the global minimum of the potential while the false vacuum is a local minimum with higher energy. We anticipate that the energy difference between these vacua produces a pressure on the domains of false vacuum causing the domain walls to collapse when this pressure becomes comparable to the surface tension of the walls [7]. Therefore, the domain wall problem could be eliminated in this scenario if domain wall networks are sufficiently short-lived that they do not survive long enough to dominate the energy density of the universe. We earlier made the self-scaling argument that domain wall energy density can be expressed as ρ dw ∝ Et −1 . Let us add an exponential suppression to this domain wall energy we have ρ dw ∝ Et −1 e −αt with corresponding density parameter, where the parameter α encodes the breaking of the symmetry and we will estimate this in the 2HDM. Again, introducing a dimensionless proportionality constant, A, specifying the number of domain walls per horizon and recalling that in write The time at which the exponential suppression dominates is determined by the parameter α after which time the domain wall density relative to critical begins to decrease. In other words, the domain wall density is maximal at t max = 1/α. Therefore, the maximum value is given by 4 The minimal theoretical requirement is that Ω max dw < 1 such that the domain walls do not dominate the energy density of the universe by the end of their scaling phase. Inserting numerical values into (4.2) we obtain the lower bound It has been established that kink solutions in the Z 2 -symmetric 2HDM have dimensionless energy of order 1 for physically viable values of physical observables [10]. Note that the time of radiation-matter equality, t eq ∼ 10 12 seconds. As such, it is apparent that even the smallest acceptable exponential suppression, α min , places the time of maximum domain wall density well within the radiation dominated epoch provided domain walls are not ultra-light, i.e. for large tan β whereÊ becomes small. It should also be noted that the collapse of these domain walls could still have undesirable effects on the Cosmic Microwave Background (CMB) and Big Bang Nucleosynthesis (BBN) [22] conflicting with current cosmological constraints on these processes. One may wish to consider constraints on domain walls arising at these epochs, however, the domain wall density (4.1) suggests that domain wall domination in the 2HDM would not have arisen by the BBN epoch. Moreover, without choosing particular combinations of physical parameters, e.g. ultra-light domain walls, such defects will have already come to dominate the universe prior to recombination. We anticipate that limits on the suppression coefficient, α, will allow constraints on the soft-breaking parameterm 2 to be obtained such that the domain walls can be made cosmologically benign. We have performed (2+1) dimensional simulations for the global scalar field theory of the 2HDM with approximate Z 2 symmetry on a regular grid of P 2 points for P = 4096 with Minkowski metric (for details of the simulation procedure see [10]). Simulations are performed in (2+1) dimensions for computational ease. Nonetheless, these simulations are a good approximation of the behaviour in (3+1) dimensions, as shown in [10]. The evolution of a set of such simulation is presented in Fig. 3 for various values of the soft breaking parameter,m 2 . In these simulations domain walls are short-lived with the entire space coming to be dominated by a single vacuum at late times. The time at which the field comes to occupy the true vacuum throughout the space is determined by the size of the symmetry breaking term,m 2 . This collapse has a significant effect on the scaling behaviour of domain walls in the approximately Z 2 symmetric 2HDM. The number of domain walls as a function of time in (2+1) dimensions are presented in Fig. 4 for various values of the dimensionless soft-breaking parameter,m 2 . The time evolution of the number of domain walls is obtained as an average over 10 realizations. Fig. 4 shows that the number of domain walls in the 2HDM with approximate Z 2 symmetry decreases much more rapidly than the t −1 scaling found in models with exact discrete symmetries [1]. The number of domain walls appears to follow an exponentially suppressed power law, N dw ∝ t −n e −αt . A non-linear least squares fitting to the number of domain walls for an exponentially suppressed power law are also included in Fig. 4 for both n = 1 and n as a fitted parameter. The exponential suppression parameter, α, exhibits the approximate relationship to the soft-breaking parameter,  Hence, in order for the exponential suppression of domain wall density to be sufficiently large to avoid domination, we obtain the limit  Assuming a small number of domain walls per horizon such that A is of order unity, this  limit suggests a small value of m 12 relative to the electroweak scale (around µeV order) would sufficiently modify domain wall scaling to avoid their domination.
In a Z 2 -symmetric Type-II 2HDM, the origin of a small effective m 2 12 parameter may be attributed to anomalous QCD instanton effects [30,31]. In the absence of a Peccei-Quinn mechanism [32,33], we may nonetheless adapt their results and conservatively estimate the size of the anomalous Z 2 -breaking parameter m 2 12 from the would-be PQ instanton potential as follows: where Λ QCD ∼ 0.3 GeV is the QCD confinement scale, n G = 3 is the number of the SM quark generations, and θ is the well-known strong CP-odd phase generated by QCD instantons. A non-zero value of θ would induce a non-zero Electric Dipole Moment (EDM) for the neutron [34]. Current experiments place an upper limit on θ < ∼ 10 −10 -10 −11 [35]. On the other hand, combining (4.8) with (4.7), we obtain a lower limit on θ, This suggests that the parameter tan β should lie in the narrow interval: 0.3 < ∼ tan β < ∼ 3.

Biased Initial Conditions
One can also avoid domain wall domination in models with an exact discrete symmetry by biasing the initial conditions such that the degenerate vacua are selected with unequal probability. In many studies of domain wall dynamics, including our own simulations of 2HDM domain walls [10], it is assumed that domain walls evolve from initial conditions where each of the degenerate vacua are selected with equal probability, i.e. that the scalar field(s) are in thermal equilibrium before the phase transition [36]. However, if the initial conditions in the early Universe have some bias towards one of the vacua, smaller domains of the disfavoured vacuum should form surrounded by larger regions where the field lies in the preferred vacuum [7]. These small domain walls should then collapse rapidly. This has been demonstrated for the Goldstone model in (2+1) and (3+1) dimensions [36]. Of course, these initial conditions must be viable for the Higgs fields in the early Universe if such biased 2HDM domain walls are to be of interest for cosmology.
Assuming an exponential suppression of domain wall scaling, the lower bound of (4.4) still holds in this case. The aim now becomes to relate this to the bias parameter, ε. We have performed (2+1) dimensional simulations with P = 4096 for the global scalar field theory of the 2HDM with exact Z 2 symmetry with Minkowski metric from biased random initial conditions. Specifically, we produce random initial conditions for the scalar fields normally distributed around ε such that one vacuum is selected with greater probability. The evolution of a set of such simulation is presented in Fig. 5. We find that domain wall are short-lived with the entire space coming to be dominated by the preferred vacuum at late times. This behaviour is qualitatively similar to that found in the softly-broken Z 2 case of Fig. 3. The number of domain walls as a function of time in (2+1) dimensions for biased initial conditions are presented in Fig. 6. The time evolution of the number of domain walls is obtained as an average over 10 realizations. Fig. 6 shows the number of domain walls decreasing with a similar profile to that seen in the case of approximate symmetry. The number of domain walls appears to follow an exponentially suppressed power law, N dw ∝ t −n e −αt . A non-linear least squares fitting to the number of domain walls for an exponentially suppressed power law are also included in Fig. 6. The exponential suppression parameter, α, shows an approximate linear relationship to the biasing, Hence, in order for the exponential suppression of domain wall density to be sufficiently large to avoid domination, we obtain the limit Again, assuming a small number of domain walls per horizon such that A is of order unity, this limit suggests a very small biasing of the initial conditions would be sufficient to avoid domain wall domination.

Conclusions
In this article we have considered the phenomenological implications of domain walls in 2HDMs with an exact or approximate Z 2 symmetry. We have obtained cosmological constraints on the Higgs masses and mixing angles such that domain walls in these models do not dominate the energy density of the Universe today. We find that domain wall domination can always be avoided for sufficiently large or small values of tan β where domain walls become ultra-light, i.e. the energy of the neutral vacuum solution tends to zero. Moreover, for Type-I 2HDMs with spontaneous breakdown of the Z 2 symmetry, we find that domain wall domination can only be avoided today for tan β > 10 5 for scalar masses larger than 100 GeV.
We have also demonstrated that domain wall networks in (2+1) dimensional simulations can be made to collapse by rendering the discrete symmetry approximate via a small symmetry breaking term. We find that the time evolution of the number of such domain walls exhibits an exponential suppression of the approximate power law scaling found in [10]. The collapse rate of the domain walls is linearly related to the soft-breaking parameter squared, m 2 12 . Consequently, we find that a soft-breaking parameter m 12 ∼ 10 −6 eV is sufficient to avoid domain wall domination by the end the scaling phase of their evolution. For a 2HDM of Type II, this suggests a corresponding lower limit on the CP-odd phase θ generated by QCD instantons, i.e. θ > ∼ 10 −11 /s β c β . This estimate is in some tension with an upper limit on θ < ∼ 10 −10 -10 −11 coming from the non-observation of a non-zero EDM for the neutron. Taking this last constraint into account, we obtain an upper and lower limit on the key parameter tan β, i.e. 0.3 < ∼ tan β < ∼ 3. We anticipate similar exponential suppression of domain wall scaling will be obtained in 2HDMs with approximate CP1 and CP2 symmetries, or in any alternative scenario which breaks the Z 2 symmetry.
Finally, we have demonstrated that domain walls evolving from biased initial conditions can be similarly short-lived with qualitatively similar behaviour to the case of an approximate symmetry. We find that domain walls in our biased simulations also experiencing an exponential suppression of their scaling. In particular, we have derived in (5.2) a lower limit on the biasing parameter ε of initial conditions such that domain wall domination can be avoided by the end of scaling. Results obtained for approximate and biased discrete symmetries can both provide means of avoiding the late-time scaling problems which domain walls ordinarily present, and hence they could be used to make 2HDMs with discrete symmetries cosmologically safe.