Simulations of domain walls in Two Higgs Doublet Models

The Two Higgs Doublet Model predicts the emergence of 3 distinct domain wall solutions arising from the breaking of 3 accidental global symmetries, Z2, CP1 and CP2, at the electroweak scale for specific choices of the model parameters. We present numerical kink solutions to the field equations in all three cases along with dynamical simulations of the models in (2+1) and (3+1) dimensions. For each kink solution we define an associated topological current. In all three cases simulations produce a network of domain walls which deviates from power law scaling in Minkowski and FRW simulations. This deviation is attributed to a winding of the electroweak group parameters around the domain walls in our simulations. We observe a local violation of the neutral vacuum condition on the domain walls in our simulations. This violation is attributed to relative electroweak transformations across the domain walls which is a general feature emerging from random initial conditions.


Introduction
Although the Standard Model (SM) of particle physics has been conspicuously successful when confronted with experimental data, it is still not satisfactory as a complete or fundamental theory of nature. Several observed phenomena in cosmology cannot be accounted for in the SM, including dark energy [1], dark matter [2] and the observed baryon asymmetry of the Universe [3,4]. This necessitates some new physics beyond the Standard Model (BSM) to explain these phenomena. There are also motivations from particle physics to extend the SM. One interesting idea was that the electroweak and strong forces are unified into some grand unified theory (GUT) at high energy [5,6]. Other BSM theories aim to address the so-called gauge hierarchy problem of the strength of forces in the SM and gravity [7], or provide mechanisms by which neutrinos can acquire mass [5,[8][9][10][11] as first 2 The Two Higgs Doublet Model

Introduction
The Higgs mechanism for electroweak symmetry breaking was confirmed at the LHC by the discovery of a Higgs boson of mass 125 GeV [27]. Properties of this scalar so far match those of the SM [28], nonetheless, current data does not rule out the existence of more scalar particles. Indeed the suggestion of a more complicated scalar sector far predates the SM Higgs discovery [14]. There have been numerous studies, both theoretical and experimental, made into the possibility of extensions to the SM involving additional scalar particles. One of the simplest models is the 2HDM where one additional complex scalar doublet is added to the SM. 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 ± [16]. The other three scalar degrees of freedom correspond to Goldstone bosons which are absorbed JHEP01(2021)105 into the longitudinal components of the electroweak gauge bosons, W ± and Z 0 . As stated earlier, there are currently several observed phenomena in cosmology which indicate the need for BSM physics including dark matter and the need for additional sources of CP violation. The 2HDM has been used extensively in the literature to explore, for example, new particle physics phenomenology beyond the Standard Model which could explain these observations [15,16,24,[45][46][47].
The Lagrangian of the electroweak-Higgs sector of the 2HDM is given by where Φ 1 and Φ 2 are the complex scalar Higgs doublets while (2.2) and are the SU(2) L and U(1) Y field strengths, respectively. The gauge-covariant derivative is given by where σ a are the Pauli matrices. The equations of motion for the Higgs and electroweak gauge fields can be derived to be (2.6) and where summation over repeated indices is implied. The right-hand sides (r.h.s.) of (2.6) and (2.7) are U(1) Y and SU(2) L currents, respectively, which are defined as and The most general form of the tree-level 2HDM scalar potential is [14,15] Table 1. Parameters in the 2HDM for which the model possesses an accidental discrete symmetry. Dashes signify the absence of a restriction on the corresponding parameter. See [30,51] for a complete classification of all possible symmetries.
Symmetry µ 2 Table 2. Parameters for the 2HDM in the reduced basis, Im m 2 12 = Im (λ 5 ) = Im (λ 6 ) = 0 and λ 6 = λ 7 , for which the model possesses one of the three accidental discrete symmetries. Dashes signify the absence of a restriction on the corresponding parameter.
The general 2HDM potential has 14 free parameters where µ 2 1 , µ 2 2 , λ 1 , λ 2 , λ 3 and λ 4 are real and m 2 12 , λ 5 , λ 6 and λ 7 are complex. It is common to impose some additional symmetry on the model, thereby reducing the number of free parameters. In our context we note that the 2HDM can possess 3 accidental discrete symmetries which predict domain wall solutions under the parameter reductions shown in table 1 [15,30,[48][49][50]. It is useful to note that one can further reduce the 2HDM parameter space whilst keeping all 3 accidental symmetries of interest using the reduction λ 6 = λ 7 and m 2 12 , λ 5 and λ 6 being real. In this reduced basis the 3 discrete symmetries are obtained via the constraints given in table 2. However, in this basis, the restrictions for a CP2 symmetry are simply a subset of those for Z 2 [51]. In fact, the CP2 model possesses an enhanced S 2 × Z 2 symmetry for µ 2 1 = µ 2 2 , λ 1 = λ 2 and λ 6 = λ 7 = 0 [52].

Vacua
The vacua in the 2HDM can be separated into 3 types [15]: (i) CP-preserving vacua where the vacuum expectation values (VEVs) have no relative phase between them, (2.11) (ii) CP-breaking vacua where the VEVs have a relative phase, ξ, (2. 12) and (iii) charge-breaking vacua where there is a non-zero upper component in one of the doublets, (2.13)

JHEP01(2021)105
The parameters v 1 , v 2 , v + and ξ are referred to as the vacuum manifold parameters. It can be demonstrated that charge-breaking vacua are the most general one would need by exploiting electroweak gauge freedom [15]. In other words one can produce a general parametrization of the Higgs doublets via an electroweak gauge transformation (EWGT) of the charge-breaking vacuum. Let us therefore parametrize the Higgs doublets via an EWGT of the charge-breaking vacuum, with where G a ≡ G a /v SM , and θ and G a are the would-be Goldstone bosons. We can write the U(1) Y current as The U(1) Y current decomposes into a term involving the gradient of ξ and terms involving covariant derivatives of the group parameters. Let us consider the gauge fields as pure gauges, obeying the gauge-fixing condition: W 0 = 0 and B 0 = 0. As a consequence, the SU(2) L and U(1) Y field strengths given in (2.2) and (2.3) vanish identically, i.e. W a µν = B µν = 0. Moreover, in pure gauge, the fields W µ and B µ can be removed from the electroweak currents, defined in (2.8) and (2.9), by an EWGT. Upon an appropriate EWGT, the electroweak currents then take on the simpler form, Since the left-hand sides (LHS) of (2.6) and (2.7) vanish, so should their RHSs. Consequently, the unitary gauge transformation matrix U in (2.15) and the spatial profile of the would-be Goldstone fields, θ and G a , cannot be arbitrary. Nevertheless, employing the equation of motion of the scalar doublets (2.5), one can show that the electroweak currents are divergenceless, i.e. ∂ µ J µ B = 0 and ∂ µ J a,µ W = 0. Given that the electroweak currents J µ B and J a,µ W are zero in the boundaries, it is obvious that these currents must vanish everywhere, at least in the approximation with one spatial dimension. Hence, the equations of motion for the electroweak gauge fields are automatically satisfied for these minimum energy kink solutions. Moreover, the vanishing of the field strength tensors, W a µν and B µν , implies that their contribution to the vacuum energy is zero. Hence, the equations of motion for the scalar sector (2.5)-(2.7) enforce minimization of the energy density as well.

JHEP01(2021)105
Taking into account the above considerations, the complete 2HDM Lagrangian of direct relevance to our study is given by (2.20) It is the equations of motion derived from this Lagrangian which are solved numerically in our simulations of domain walls in the next sections.

Parametrizing the global field theory
An alternative representation of the scalar potential is the so-called bilinear field-space formalism [30,48,53,54], where the first term contains the quadratic mass parameters and the second contains the quartic couplings. The specific forms of the 4-vectors and tensor are In (2.22a), Φ ≡ (Φ 1 , Φ 2 ) T and σ µ are the Pauli matrices including the identity, σ 0 = I 2 .
Note that in (2.22c) we have introduced the short-hand notations It should be noted that the 4-vector, R µ , is invariant under a unitary transformation and can, therefore, be equally written in terms of the vacuum manifold parameters: This allows us to relate the components of R µ to the vacuum manifold parameters via the expressions

JHEP01(2021)105
This general parametrization admits the possibility of charge-violating solutions. One can ensure that a solution respects charge symmetry by checking/demanding that R µ have null norm, i.e. R µ R µ = 0. This is the so-called neutral vacuum condition [48] and is achieved by having v + = 0. It is apparent that there are several ways one can choose to represent the Higgs doublets in the 2HDM. Possibly the most obvious representation is the so-called linear representation, where the components of the doublets are parametrized by 8 real scalar fields, φ i . Generating a general parametrization of the doublets via a local electroweak transformation of the general vacuum (2.13) with the U(1) Y × SU(2) L matrix, , (2.26) one produces an alternative 8-parameter representation of the doublets involving the vacuum manifold parameters, v 1 , v 2 , ξ and v + , and the electroweak group parameters γ 1 , γ 2 , γ 3 and θ. One could equally parametrize the doublets in terms of the components of R µ and the 4 electroweak parameters γ 1 , γ 2 , γ 3 , θ. We can relate the vacuum manifold parameters back to the linear representation via the expressions for R µ in terms of the scalar fields φ i . One can also introduce the SU(2) L invariant object Φ T 1 iσ 2 Φ 2 and promote R µ to a 6-vector [30,51], allowing us to exchange the group parameter θ in favour for R 4 and R 5 , Hence, we can express the remaining electroweak parameters as These parametrizations are summarized in table 3.

Simulation techniques
In what follows we will want to perform two distinct tasks. The first is to calculate the minimum energy solutions to a set of equations (either the full 8 fields, or a restricted ansatz). This is done using the Gradient Flow technique. The other is to evolve the full equations of motion including the 2nd order time derivatives.

JHEP01(2021)105
Linear Representation Table 3. Alternative parametrizations of the vacuum manifold in the 2HDM both in the original scalar field space and bilinear R-space. For definitions and relations among parameters, see text.

Gradient flow
Minimum energy field configurations for the full 8 field case, ie. φ 1 , .., φ 8 , can be obtained via variational techniques as the solution to the set of ordinary differential equations, written here for one spatial dimension, appropriate for search for kink-like solutions. This can be done by introducing an artificial simulation time, t, and solve the gradient flow equations (see, for example, [30]), The solution to these equations gives the minimum energy solution in the long-time limit subject to appropriate boundary conditions where which represent 8 coupled 1st order differential equations that can be discretized using standard techniques. We will refer to solutions generated in this way as relaxed solutions.
The gradient flow equations are evolved on a discrete grid where both spatial and temporal derivatives are approximated to second-order. We will also do this for various ansatzes for the fields, φ n . The approach is exactly the same but typically this will reduce the dimensionality of the problem, and hence the number of coupled equations.

Full dynamics in two and three dimensions
Although minimum energy kink solutions obtained via gradient flow are useful for investigating properties of domain walls such as their width and energy, more complex simulations of domain wall dynamics are required to determine features such as scaling behaviour and related domain wall interaction properties. The equations of motion for the global scalar field theory are discretized and solved on a regular grid in (2+1) and (3+1) dimensions. Simulations begin from normally distributed random initial conditions in each field of the linear representation, φ 1 , .., φ 8 . These initial conditions are chosen to broadly represent the outcome of a 2nd order phase transition by randomly assigning the field configuration to one of the degenerate vacua at each grid point. These initial conditions are not intended JHEP01(2021)105 to model the detailed dynamics of realistic phase transition, rather, they are a minimal practical choice. They initially produce large energy gradients between grid points. We smooth these discontinuities over adjacent grid points by introducing an artificial damping term (φ) with coefficient, ε [38][39][40]. This damping term is turned off after an appropriate, small number of time steps to restore the physical equations of motion. This numerical technique produces some transient effects at the point when the equations transition to the physical regime, however, this is short-lived. In results presented later, we choose ε = 0.5 for the first 200 time steps before the damping effect is removed unless stated otherwise. Spatial derivatives are approximated to fourth order and temporal derivatives to second order. Our simulations use periodic boundary conditions. Therefore, we must make sure to only consider results up to the so-called light-crossing time, τ = 1 2 P ∆x, where P is the number of spatial grid points and ∆x is the grid spacing. Beyond the light-crossing time the domain walls could have crossed the periodic boundaries and interacted with themselves making any further dynamics unphysical. Note that we have applied the numerical techniques described here to the Goldstone model and obtained the well-established ∝ t −1 scaling that one would expect.
We also investigate the evolution of domain wall networks in an expanding universe. However, the thickness of a domain wall is constant for a given set of model parameters, whereas the grid spacing will increase with time. Therefore, in an expanding universe their comoving thickness will decrease in inverse proportion to the scale factor, a, [55]. This results in resolution issues as the lattice expands, with the width of the kink becoming smaller than the grid spacing. In order to ameliorate this effect, one modifies the equations of motion to give a constant comoving resolution, where α and β are constants. 1 This is the so-called PRS algorithm [44]. Setting β = 0 produces a fixed comoving thickness for the domain walls and α+β/2 = N yields unchanged dynamics for domain walls in (N + 1) dimensions [44,55]. In the results presented later, simulations start from an initial conformal time η 0 = 1 and use a power law, a = η γ , for the scale factor. This algorithm has been widely used in studies of domain walls (see, for example, [41,56]). In the literature, the PRS algorithm [44] has mostly been applied to simple scalar field theories such as the Goldstone model. Even though the PRS algorithm has not been tested before for more elaborate models, we will nevertheless assume its applicability to the 2HDM that we will be studying here. The validity of such a consideration is firmly supported by the fact that our results presented in the following sections will exhibit a qualitatively similar behaviour to that obtained in simulations using the Minkowski metric where the PRS algorithm is absent.

Z 2 symmetry
A Z 2 transformation of the Higgs doublets is given by (4.1) The Z 2 -symmetric 2HDM potential is obtained by applying the parameter restrictions of table 1 on the general 2HDM potential, (2.10), and can be written as For the CP-preserving scenario given in (2.11), its VEVs can be calculated in terms of the potential parameters [30], In order to have a stable vacuum there must be an all-positive Higgs spectrum, i.e. we must have real, non-negative eigenvalues for all 5 physical Higgs states. It is most convenient to compute the mass matrices using the representation, where ϕ + i are complex scalar fields and the CP-even/odd fields can be expressed as respectively, where c α = cos α and s α = sin α, and similarly for c β and s β . It should be obvious that the representation (4.4) reduces to (2.12) in the vacuum. Since the Z 2symmetric 2HDM is CP conserving, there is no mixing between the CP-even and CP-odd states and ξ = 0. Therefore, the mass matrix for the neutral sector becomes block diagonal. The mass matrices in the Z 2 -symmetric 2HDM are found to be (4.6) (4.8) From the above matrix relations, it is not difficult to calculate the physical squared masses for all scalar fields, i.e.
along with 2 zero eigenvalues corresponding to the would-be Goldstone bosons, G 0 and G ± .
One can re-express the scalar potential in terms of these masses and the CP-even/odd mixing angles α/β, where we have also rescaled for dimensionless energy per unit area.

Kinks solutions
In [30] kink solutions were obtained as a function of the vacuum manifold parameters with (2.11) and (2.12) assumed as an ansatz reducing the number of gradient flow equations to two. Despite this ansatz being well-motivated, allowing one to obtain minimum energy solutions, it only allows exploration of a limited region of the field configuration space. Hence, the solutions available are restricted from the outset. In contrast, we make no assumption of the vacuum manifold parameterization in this work and obtain relaxed solutions from the general field configuration and the general vacuum ansatz (2.13). The general energy density of the Z 2 -symmetric 2HDM in one dimension is given by In the linear representation, the gradient flow equations which minimize the energy per unit area are given by

JHEP01(2021)105
The energy density in terms of the vacuum manifold parameters is, (4.16) The gradient flow equations for the vacuum manifold parameters are then given by Note that the relations between the 2HDM potential parameters and physical parameters can be found in appendix A.
The relaxed solutions to (4.15) and (4.17) are presented in figure 1. It can be seen that the solutions to these two systems of equations are equivalent, having the same energy and identical profiles for the components of R µ , for specific model parameters, and this is the case for all parameter choices which we have tried. Moreover, the relaxed solution in the full field configuration of figure 1 can be brought to the same form as the solution for the vacuum manifold parameters via a global transformation of the fields multiplying by U ∈ SU(2) L × U(1) Y -a global electroweak transformation. In other words, there are a myriad of equivalent minimum energy solutions in the general parameterization which are related via global electroweak transformations. One of these is the solution found in [30] which used the restricted ansatz with ξ ≡ 0 and v + ≡ 0, and indeed the minimum energy solutions of (4.17) have this property.
In both cases a kink forms in R 1 and both have the same vacuum values for the components of R µ given by (2.23). The quantities R 1 and R 2 are odd under a Z 2 transformation whilst R 0 and R 3 are even. We may then define a gauge-invariant topological current, and corresponding topological charge, Note that Q = 1 corresponds to a kink and Q = −1 to an anti-kink.

Dynamical simulations from random initial conditions
We have performed (2+1) dimensional simulations with P = 1024 and P = 4096 up to t = τ for the global field theory of the Z 2 -symmetric 2HDM with Minkowksi metric, and JHEP01(2021)105 We find that there is a local violation of the neutral vacuum condition, R µ R µ = 0, on the domain walls in our simulations as shown in figure 2; a feature of the field dynamics not found in the minimum energy solutions, indicating that there is something more complicated going on than just evolution of a network of minimum energy configurations. We will return to this observation later in our discussion. Equivalent results for a domain wall network in a radiation dominated FRW universe are presented in figure 4. Again, we find a local violation of the neutral vacuum condition in the vicinity of the domain walls.
The aforementioned dynamical features are manifest for a wide range of physical parameters. To illustrate this we consider a benchmark scenario motivated by the socalled Maximally Symmetric 2HDM [50,57], which has a common heavy Higgs mass, [50,57], we assume here a Type-I Yukawa realization, so that the stringent constraints, mainly due to B-meson observables, can be avoided [58,59]. For such a scenario, the domains and spatial variation of R µ R µ are presented in figure 5.
The number of domain walls, N dw , as a function of time in (2+1) dimensions with Minkowski and FRW metrics in radiation and matter dominated eras are presented in figure 6. The time evolution of the number of domain walls in the system is obtained as an average over 10 realizations. In both Minkowski space and an FRW radiation era we find domain walls in the Z 2 -symmetric 2HDM do not scale in the standard way, whereas in the FRW matter era it is compatible with N dw ∝ t −1 . There is an obvious trend that the faster the expansion rate, the less the domain wall network deviates from scaling. It should be noted that the matter era simulations are beginning to deviate and we anticipate the deviation would become more pronounced in simulations of greater dynamic range. Figure 6 JHEP01(2021)105  shows a deviation from the t −1 scaling found in simpler models such as the Goldstone model [34]. In the 2HDM we find more domain walls at late times than one would naively expect. This can be seen visually in that the number of walls from panel to panel in all of figures 2, 4 and 5 clearly does not decrease by a factor of two as one progresses from left to right. This non-standard scaling is akin to the effect found in the so-called kinky vorton model [39]. In the kinky vorton model a Noether current produces a conserved charge which condenses on the domain walls. The winding of the condensate field around the walls slows the collapse of the domain walls modifying their scaling behaviour. Our results suggest that there may be some extra interactions amongst the fields which prevents the walls from decaying as quickly as causality permits.   the Z 2 -symmetric 2HDM is given in figure 7. As in the case of (2+1) dimensions, we observe a local violation of the neutral vacuum condition in the vicinity of the domain walls as shown in figure 8. We also find a deviation from the expected power law scaling as shown in figure 9, albeit with significantly less dynamic range compared to the case of (2+1) dimensions.

Kink interactions
The deviation from the expected power law scaling observed in figures 6 and 9 may suggest some interaction among the scalar fields which prevents the domain walls from scaling as fast as they otherwise would. This could be attributed to the fact that the kink-antikink interaction can be attractive or repulsive depending where they form in the field configuration. To investigate this we calculate the force between a kink and an anti-kink in (1+1) dimensions with a relative SU(2) L × U(1) Y global rotation between them. The energy-momentum tensor for the 2HDM is given by Therefore, in (1+1) dimensions, where dots and primes denote temporal and spatial derivatives, respectively. The momentum of the field configuration over some semi-infinite length is then (4.22) and the corresponding force Using the equations of motion to substitute forΦ i andΦ † i in (4.23) and integrating, we find (4.24) In order to evaluate the lower limit, x = b, we can derive asymptotic expressions for the kink and condensate. (with i = 1, 2), the Euler-Lagrange equations reduce to We consider the force on a field configuration composed of a kink-anti kink pair with a relative SU(2) L × U(1) Y transformation between the kink and anti kink as follows whereΦ 1,2 (x) are static solutions 3 of the type seen in figure 1, Φ 0 1,2 are the CP-preserving vacua of (2.11) and is a global SU(2) L × U(1) Y matrix. It should be noted that this is an approximate solution, becoming exact in the limit that U becomes an identity matrix. Nonetheless, this field configuration has the desired profile over the range of integration. For −a b a, we JHEP01(2021)105 haveΦ i (x + a) − Φ 0 i 0 and, hence,Φ i (x + a) 0, allowing us to linearize (4.24) in these small quantities. Firstly, we can linearize the spatial derivatives in (4.24) as where we have defined the short-hand notationΦ i,± =Φ i (x ± a). The 2HDM potential of (4.2) linearizes as Hence, the linearized version of the force on the field configuration can be written as, (4.31) Since the static solutions are asymptotically flat, their gradients vanish at infinity. Therefore, there is no contribution from the upper limit of (4.31). Furthermore, substituting the explicit form of the doublets, we reduce the above expression to In order to evaluate the lower limit, let us define the asymptotic expressions where µ = √ 2µ 1 , ν = √ 2µ 2 and ρ and σ are normalization factors that are undetermined at this level of approximation. Using these asymptotic expressions for the kink solutions we can now write the force as where R = 2a is the kink separation. Therefore, the interaction energy is We see that the group parameters, A and θ, are present in the expression. This indicates that the sign of the force will change depending on the values of the group parameters. As such, there can be repulsive as well as attractive interactions between domain walls in different parts of the field configuration. Note that for U = I (where θ = 0 and A = 1) we recover the usual negative interaction energy signifying an attractive force between the kink and anti kink. Now imagine an approximately circular wall evolving as part of the wall network in (2+1) dimensions (or an approximately spherical wall in (3+1) dimensions). The standard picture would be for the walls to collapse as fast as causality would allow under their own tension -this is what leads to the expectation that N dw ∝ t −1 . However, in this case we have shown that there are extra forces associated with the relative phases of the vacuum which could be, and indeed will be different on each part of the wall. These forces will clearly interfere with the dynamics of the wall network and we can expect that the scaling behaviour will be modified.

Relatively gauge-rotated vacua
Let us now consider the effect of introducing a relative electroweak transformation between vacua on either side of a kink obtained via gradient flow. This is something which is not possible in the case of the standard domain wall solution in the discrete Goldstone model, but which is manifest in the 2HDMs with discrete symmetries, and would be expected in simulations from random initial conditions. In order to achieve this, we introduce a relative transformation between the vacuum field configuration at the boundaries and impose Dirichlet boundary conditions on the gradient flow, where U has the form of (2.26). If we consider the effect of each group parameters separately, we see that a non-zero γ 1 excites the upper components of the CP-preserving vacua, while θ and γ 2 excite the imaginary part of the lower components.
In order to investigate the impact of U = I, we have first calculated solutions for a range of the parameter γ 1 keeping γ 2 = γ 3 = θ = 0 fixed throughout. The peak value of R µ R µ as a function of γ 1 is presented in figure 10 along with the corresponding energy per unit area of the solutions. The solution with lowest energy is indeed the solution where there is no relative transformation between the boundaries. For this solution we observe no violation of the neutral vacuum condition as expected.
However, for any other value of γ 1 we obtain higher-energy solutions where R µ R µ is non-zero and peaked locally at the position the kink as seen in our simulations from random initial conditions. An example solution for γ 1 = π/2 and the corresponding R-space profiles are presented in figure 11. We note that the neutral vacuum condition is respected only when the upper components of the doublets are zero at the rotated boundary, i.e. in the absence of v + . Of course, the violation of neutral vacuum and the increased energy of such solutions as shown in figure 10 are not unrelated. Since the solutions which satisfy the minimum energy condition are those which have v + = 0, in order to remove the additional energy in the neutral vacuum-violating solutions we would need to eliminate the locally non-zero v + . In this case we obtain the neutral vacuum solution. We also observe an increase in the energy per unit area as in figure 10 for solutions with transformations involving the other group parameters, γ 2 , γ 3 and θ.
We have also investigated solutions when θ and γ 2 = 0. Since they both excite the same components of the field configuration (φ 4 and φ 8 ) with opposite sign, we only present results for an example solution for θ = π/2 with γ 1 = γ 2 = γ 3 = 0 which is presented in figure 12.   Figure 11. Kink solution with γ 1 = π/2 at the right hand boundary (left) and corresponding R-space profiles for M H = M A = M H ± = 200GeV, tan β = 0.85 and cos(α − β) = 1.0. Note that there is no obvious wall solution in any of the fields φ 1 , .., φ 8 , but there is clearly a kink in the profile of R 1 . Moreover, in this case R µ R µ is non-zero at the position of the kink, and R 4 is also non-zero, but negative. Since θ = 0, R 4 is non-zero while R 5 = 0 (see (2.27)) and we also have that around the kink not observed in the minimum energy solution of figure 1. Specifically, a non-zero γ 1 leads to a solutions with non-zero R 4 and R 5 on the kink and hence a violation of the neutral vacuum condition. For non-zero θ (or γ 2 ) we obtain solutions with a non-zero R 2 around the kink. Corresponding plots of the vacuum manifold parameters are given in figure 13. For non-zero γ 1 we obtain a non-zero v + peaked at the kink indicating a local violation of the neutral vacuum condition. For non-zero θ we observe a non-zero ξ around the kink. Spatial variation of the vacuum manifold parameters and group parameters from a (2+1) dimensional simulation are presented in figures 14 and 15, respectively. The profiles of the vacuum manifold parameters are consistent with those obtained for the kinks with rotated vacua. Specifically, the peak in v + and the profile of ξ seen in figure 13 are manifest in the simulation results. We also observe a winding of the parameters θ, γ 2 and γ 3 around the domain walls. If these features are indeed akin to the kinky vorton model, one would expect there to be some Noether-like charge winding around the domain walls.
A plot showing the time evolution of θ is presented in figure 16. We see that there is some clockwise winding of the electroweak angle, θ, around the domain wall. Note that figure 16 shows a select spatial region containing a circular domain wall at late times in our simulations.

CP1 symmetry
A CP1 transformation 4 of the Higgs doublets is given by Using the parameter restrictions of table 1, the CP-symmetric 2HDM potential can be written as For the CP-breaking vacuum (2.12), the VEVs are calculated in terms of the potential parameters in [30]. The final results are given by

JHEP01(2021)105
and ξ = arccos In order to find the masses of the physical states, we must again calculate the eigenvalues of the Hessian matrix. We use the minimization conditions, and in order to further simplify the mass matrices given in [16]. In this way, the charged sector mass matrix simplifies to As expected, the charged mass matrix (5.9) has one zero eigenvalue corresponding to the Goldstone bosons, G ± , and the mass of the charged Higgs bosons is Unlike in the Z 2 case, the CP1-symmetric 2HDM has a non-zero phase, ξ, and the mass matrices for the CP-even and CP-odd scalars do not decouple as the model is not explicitly CP-conserving. Working in the {ϕ 1 , ϕ 2 , A, G 0 } basis, the G 0 axis decouples from the other 3 and, as one would expect, has a zero eigenvalue. Hence, we obtain a symmetric 3 × 3 neutral mass matrix, M 2 N , with elements, The masses for the three neutral Higgs bosons, M H 1,2,3 , are obtained as roots of a cubic characteristic equation, exactly as done in [16]. Which mass eigenstate (i = 1, 2, 3) corresponds to which Higgs particle is then determined by its relative coupling to the Z 0 boson, where O is an orthogonal 3×3-dimensional matrix that diagonalizes the scalar mass matrix M 2 N given in (5.11). In particular, the SM-like scalar will be that which has g H i ZZ close to 1.

JHEP01(2021)105
Ideally, to simulate the CP1-symmetric 2HDM we would obtain expressions for the scalar masses, M H 1,2,3 , of (5.11) and invert this system of equations to obtain expressions for the parameters of the scalar potential in terms of the physical observables as was performed for the Z 2 -symmetric 2HDM in appendix A. While this calculation is achievable in the Z 2 case, the CP1-symmetric 2HDM potential is significantly more complicated due to the mixing of CP-even/odd states requiring a more complicated system of mixing angles [60] to diagonalize (5.11). Nonetheless, the values of the neutral scalar mass eigenvalues can be obtained numerically for a given set of potential parameters. To achieve this we first simplify (5.11) by working in the reduced basis of table 2 with additional restrictions λ 4 = λ 6 = 0. Under these additional restrictions the quartic coupling parameter λ 5 fixes the charged Higgs mass as seen in (5.10). To obtain maximal spontaneous CP violation, we specify ξ = π/4 and choose tan β with v 1 = c β v SM and v 2 = s β v SM allowing us to infer the quadratic parameters, µ 2 1 , µ 2 2 and m 2 12 , from the minimization conditions (5.6)-(5.8). We choose λ 3 such that the mass matrix element (M 2 N ) 12 of (5.11) is small to achieve SM alignment. The remaining parameters λ 1 and λ 2 are used to fix the scalar masses, M H 1,2,3 .

Kinks solutions
As for the Z 2 -symmetric 2HDM, we present a generalized treatment of kinks in the CP1symmetric 2HDM from [30]. We obtain relaxed solutions from the general field configuration and general vacuum ansatz (2.13).
Inserting the general vacuum, (2.13), into the potential, (5.2), one finds the following expression for the energy per unit area for the CP1-symmetric 2HDM: The resulting gradient flow equations are These gradient flow equations simplify in the reduced basis with λ 4 = λ 6 = 0 for which we obtain physical parameters. The relaxed solutions of this system of equations in the reduced basis are given in figure 17 in both the linear representation of the Higgs doublets and for the vacuum manifold parameters, v 1 , v 2 , v + and ξ. It can be seen that a kink forms in R 2 and correspondingly in ξ. These two systems of equations are equivalent having the same energy per unit area and profiles for R µ . This is the case for all model parameters we have considered including arbitrary parameters sets in the non-reduced basis of table 1. Therefore, as in the Z 2 case, there are a myriad of equivalent field configurations related by global electroweak transformations. One such solution was obtained in [30] using the restricted ansatz v + ≡ 0. The minimum energy solution to (5.14) indeed has this property. and R 3 are even. Hence, we may define the gauge-invariant topological current,

JHEP01(2021)105
and corresponding topological charge, Again, Q = 1 corresponds to a kink and Q = −1 to an anti kink.

Dynamical simulations from random initial conditions
We have performed (2+1) dimensional simulations with P = 1024 and 4096 up to t = τ for the global field theory of the CP1-symmetric 2HDM with Minkowski, and FRW metric in the radiation and matter dominated eras using the PRS algorithm. The evolution of a Minkowski simulation is presented in figure 18 along with the corresponding spatial variation of v + , which would be zero in the case of a neutral vacuum solution, v + ≡ 0.
The domains are mapped according to the sign of R 2 as one would expect given the kink solutions in figure 17. As in the Z 2 case, we find a local violation of the neutral vacuum condition, R µ R µ = 0, on the domain walls in our simulations as shown in figure 18. The magnitude of this violation is small for physical parameter sets we have considered. Nonetheless, the violation is present and is qualitatively similar to that found in the Z 2 -symmetric 2HDM. For clarity of presentation we have also performed simulations in the reduced basis for an arbitrary set of scalar potential parameters as presented in figure 19. In this case the violation of the neutral vacuum on the domain walls which emerges in our simulations can be more easily seen. For comparison to figure 18, the corresponding physical observables for the parameters chosen in figure 19 are M 1 = 123 GeV, M 2 = 317 GeV, M 3 = 354 GeV, M H ± = 123 GeV, ξ = 1.48 rad and tan β = 1.0. It should be reiterated that the qualitative features presented here are general to all parameter sets we have considered.
The number of domain walls as a function of time in (2+1) dimensions with Minkowski and FRW metrics in radiation and matter dominated eras are presented in figure 20. The time evolution of the number of domain walls is again obtained as an average over 10 realizations. As in the Z 2 -symmetric case, we find that in both Minkowski space and FRW radiation era domain walls do not follow standard scaling, whereas the FRW matter era is compatible with N dw ∝ t −1 .

JHEP01(2021)105
Making use of these expressions, the elements of the charged mass matrix can be written as which has one zero eigenvalue for the Goldstone bosons, G ± , and a charged Higgs mass, Moreover, the elements of the 3 × 3 squared mass matrix read As in the CP1 case, the three squared mass eigenvalues of the above matrix, M 2 h,H and M 2 A , are obtained by solving a cubic characteristic equation. However, in the diagonally reduced basis of table 2, all scalar masses of the CP2-symmetric 2HDM take on a simple form given by Thus, following a rescaling for dimensionless energy, we can write the CP2-symmetric 2HDM in the diagonally reduced basis in terms of the Higgs masses, Notice that in the diagonally reduced basis, it is easier to see that CP2 becomes a restricted scenario of the Z 2 symmetry with tan β = 1 and g hZZ = 1 [52]. Therefore, without loss of generality, in all our numerical results presented in the following sections, we perform our simulations in the diagonally reduced parameters space of the CP2-symmetric 2HDM.

Kinks solutions
As for the Z 2 and CP1 cases, we have performed a generalized treatment of kinks in the diagonally reduced CP2-symmetric 2HDM. We obtain relaxed solutions from the general field configuration and general vacuum ansatz, (2.13).
Inserting the general vacuum, (2.13), into the potential, (6.2), one finds the energy per unit area for the diagonally reduced CP2-symmetric 2HDM: The resulting gradient flow equations are (6.13) In figure 21, we present the relaxed solutions both in the linear and bilinear field-space representations. In the latter representation, we see the formation of a kink along the R 1 direction. Under a CP2 transformation, the quantities R 1,2,3 are odd in a general weak basis, whilst R 0 is even. To distinguish CP2 from the other cases, Z 2 and CP1, one could in principle define the gauge-invariant topological current by means of R 3 , i.e. 14) and the respective topological charge, As before, we stipulate that Q = 1 describes a kink and Q = −1 an anti-kink. However, as we will see in the next section, in the diagonally reduced basis, we have that R 2,3 (±∞) = 0, since the permutation symmetry S 2 , under the exchange Φ 1 ↔ Φ 2 , remains unbroken. As a consequence, the definitions for the topological current and charge in (6.14) and (6.15) become singular. For this reason, we will use the definitions (4.18) and (4.19) from the Z 2 -symmetric case to identify kink and anti-kink solutions.

Dynamical simulations from random initial conditions
We have performed (2+1) dimensional simulations with P = 1024 and P = 4096 up to t = τ for the global field theory of the diagonally reduced CP2-symmetric 2HDM with JHEP01(2021)105 Minkowski metric, and FRW metric in the radiation and matter dominated eras using the PRS algorithm. The evolution of a Minkowski simulation is presented in figure 22 along with the corresponding spatial variation of R µ R µ . The domains are mapped using to the sign of R 1 given the solution in figure 21. As in the Z 2 and CP1 cases, we find a local violation of the neutral vacuum condition, R µ R µ = 0, on the domain walls in our simulations of the diagonlly reduced CP2-symmetric 2HDM as shown in figure 22.
The number of domain walls as a function of time in (2+1) dimensions with Minkowski and FRW metrics in radiation and matter dominated eras are presented in figure 23. The time evolution of the number of domain walls is again obtained as an average over 10 realizations. Both in Minkowski space and in FRW radiation-and matter-dominated eras, we find a similar feature that domain walls do not follow standard scaling in the diagonally reduced CP2-symmetric 2HDM.

Conclusions
The 2HDM can in general predict a variety of topological defects when accidental symmetries are broken. Specifically, for SU(2) L × U(1) Y preserving symmetries, there are three domain wall solutions, two vortex solutions and one global monopole solution [30]. Here, our attention has been focused solely on 2HDMs that produce domain walls after the spontaneous breakdown of the discrete Z 2 , CP1 and CP2 symmetries. Firstly, we have presented minimum energy kink solutions obtained via gradient flow. We have shown that these kinks can be visualized in the components of R µ regardless of the specific field configuration in the linear representation. For each case, we have been able to identify a gauge-invariant and conserved topological current, and the corresponding topological charge associated with these kink solutions.
For the 2HDMs with accidental Z 2 , CP1 and CP2 symmetries, we have performed numerical simulations both in (2+1) and (3+1) dimensions. In all three cases, simulations produce a network of domain walls which may manifest themselves in the components of the gauge-invariant R µ -vector, and are in agreement with changes of the associated topological charges. Specifically, domain walls in the Z 2 -symmetric 2HDM are found to flip the sign of R 1 , domain walls in the CP1 case flip the sign in R 1,2 and the CP2 case has domain walls where all R 1,2,3 change sign.
In all three domain wall-forming 2HDMs under investigation for both Minkowski and radiation-dominated FRW metrics, we have observed a deviation from the ∝ t −1 power JHEP01(2021)105 law scaling found in simpler models. In all three cases, the matter era simulations exhibit qualitatively similar profiles for the number of domain walls. We anticipate the deviation in these models will become more pronounced at later times which could be probed by simulations with greater dynamic range (cf. profile of figure 9 for (3+1) dimensions where dynamical range is reduced). Moreover, for CP2 we do indeed see inconsistency with scaling for all 3 cases in figure 23 where the deviation is more pronounced. The observation of this deviation from scaling in both Minkowski and FRW simulations shows that this deviation cannot simply be a numerical artefact of the PRS algorithm. For all three discrete symmetries, we find more domain walls at late times in our simulations than one expects. This feature can be attributed to the winding of the electroweak group parameters around the walls slowing their collapse. This feature is akin to that seen in the kinky vorton model [39]. We also observe a local violation of the neutral vacuum condition in our numerical simulations which we attribute to the relative EW rotation between domain walls. This feature was further investigated by minimizing the energy per unit area of a field configuration with a relative EW rotation between the boundaries. This minimization yields solutions with a local violation of the neutral vacuum conditions as seen in the random simulations.

JHEP01(2021)105
It would be interesting to extend the present work and discuss constraints on the theoretical parameters that are obtained exclusively from the non-observation of domain walls today. Therefore, in an upcoming paper, we plan to systematically address the particle physics phenomenological implications for all three 2HDMs with spontaneously broken discrete symmetries that result from the absence of domain walls in the present epoch.