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, $Z_2$, 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 established by the measurement of neutrino oscillations [12,13]. Therefore, it is important to consider what extensions of the SM can provide compelling explanations for these observations whilst remaining consistent with current experimental measurements.
A minimal and theoretically well-motivated extension of the SM is the so-called Two Higgs Doublet Model (2HDM) first suggested in 1973 [14], where the scalar sector of the SM is extended to include one additional complex scalar doublet (for a review, see [15]). The 2HDM has been studied extensively in the contexts of CP violation [16][17][18][19], dark matter and electroweak baryogenesis [20][21][22][23][24]. The Higgs mechanism for electroweak symmetry breaking was verified at the LHC in 2012 by the observation of a Higgs boson [25,26] with a mass of 125 GeV [27]. Properties of this scalar sector so far match those of the SM [28,29]. However, current experimental measurements do not rule out the existence of additional scalar particles. The 2HDM predicts five physical scalar particles: two neutral CP-even states, h and H, one of which may be identified with the SM Higgs; a CP-odd neutral state, A, and two charged states, H ± .
Furthermore, it is already well known [30] 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 possess under certain parameter choices. Symmetries which are currently broken, such as the electroweak symmetry, should be restored at sufficiently high energies [31] when the system exceeds the relevant critical temperature [32]. These broken symmetries are no longer directly observable. However, they should have been intact in the early Universe when temperatures were far higher than at present [31]. This suggests that the Universe has undergone a series of symmetry breaking phase transitions in its history. These phase transitions can leave relic topological defects which could be used to probe high energy physics in the early Universe [31,33,34].
Our focus in this paper will be on discrete symmetries of the 2HDM, whose spontaneous breaking predicts domain wall solutions. In particular, we are interested in numerical simulations of such topological defects by solving the relevant equations of motion [31,[35][36][37]. The results of these simulations can guide cosmological searches and reveal cosmological implications for new models providing phenomenological constraints on BSM theories complementary to those coming from particle physics experiments.
Domain walls are topological defects which emerge due to the breaking of discrete symmetries [31]. The breaking of a discrete symmetry results in a vacuum manifold containing topologically disconnected points corresponding to degenerate vacua. Phase transitions producing domain walls occur at a finite rate, therefore, causally disconnected regions of space can select different vacua. This divides the universe into so-called domains the interfaces between which are domain walls [32]. Defects which emerge from the breaking of a global symmetry such as domain walls enter a regime of dynamical scaling such that the number of defects is constant per Hubble horizon [34]. Domain walls follow a power law scaling with an exponent close to -1 as has been shown in simulations of the Goldstone model [31,36]. This is due to the fact that the walls have a surface tension under which they collapse as quickly as causality permits.
Domain walls would present an undesirable fate for the Universe were they to exist in nature. The energy density of matter and radiation in their respective epochs of domination decrease proportionally to (time) −2 . However, the energy density of domain walls decreases proportionally to (time) −1 [38]. This means that the energy density of domain walls grows relative to matter and radiation in their respective epochs and, therefore, dominates the energy density of the Universe at some late time [39][40][41]. This is the so-called domain wall problem. If cosmic domain walls are to be realized in nature, constraints are placed on domain wall-forming models such that domain wall domination does not occur [38] or, at least, occurs after present day.
The Two Higgs Doublet Model predicts three domain wall solutions which arise from the breaking of accidental Z 2 , CP1 and CP2 symmetries. We have obtained 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 all three symmetries, we find a stable network of domain walls which deviates from power law scaling in both Minkowski and Friedmann-Robertson-Walker (FRW) simulations. Most remarkably, we observe a local violation of the neutral vacuum condition on the domain walls in our simulations.
The remainder of this article exhibits the following structure. In Section 2, we outline the scalar sector of the 2HDM, derive the relevant equations of motion and discuss the conserved electroweak currents. In addition, we discuss the different types of electroweak vacua, as well as possible accidental discrete symmetries which lead to domain wall solutions. In Section 3, we give a short overview of our simulation techniques used to analyze the evolution of domain-wall networks resulting from the breaking of the accidental symmetries, Z 2 , CP1 and CP2. Thus, Section 4 presents results for Z 2 domain-wall simulations starting from random initial conditions. As explicitly demonstrated in Appendix A, our analysis was aided by a suitable parametrization of the 2HDM potential in terms of the physical scalar masses and their mixings. Likewise, Sections 5 and 6 present kink solutions and results of domain wall simulations from the breaking of CP1 and CP2 accidental symmetries, respectively. For all three possible discrete symmetries, Z 2 , CP1 and CP2, we discuss the dynamical features of the emerging domain walls, including their cosmological scaling. Finally, Section 7 summarizes the main results of our study and presents an outlook for future work.

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 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,[42][43][44].
The Lagrangian of the electroweak-Higgs sector of the 2HDM is given by where Φ 1 and Φ 2 are the complex scalar Higgs doublets while 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 where summation over repeated indices is implied. The right-hand sides (RHS) of (2.6) and (2.7) are U(1) Y and SU(2) L currents, respectively, which are defined as 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,48] for a complete classification of all possible symmetries. 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,[45][46][47]. 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 [48]. In fact, the CP2 model possesses an enhanced S 2 × Z 2 symmetry for µ 2 1 = µ 2 2 , λ 1 = λ 2 and λ 6 = λ 7 = 0 [49].

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, (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, 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, 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.
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,45,50,51], 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 λ ab = λ a + λ b ,λ ab = λ a − λ b , λ abc = λ a + λ b + λ c andλ abc = λ a + λ b − λ c . 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 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 [45] 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, 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,48], 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.
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 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, ε [35][36][37]. This damping term is turned off after an appropriate, small number of time steps to restore the physical equations of motion. 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.
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, [52]. 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 [41]. Setting β = 0 produces a fixed comoving thickness for the domain walls and α+β/2 = N yields unchanged dynamics for domain walls in (N + 1) dimensions [41,52]. In the results presented later, simulations start from an initial conformal time η 0 = 1 and use a power law, a = η γ , for the scale factor.

Z 2 Symmetry
A Z 2 transformation of the Higgs doublets is given by 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], 1 Note that a dot, here, signifies differentiation with respect to conformal time.
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 α/β,

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 (4.14) In the linear representation, the gradient flow equations which minimize the energy per unit area are given by The energy density in terms of the vacuum manifold parameters is,  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 Fig. 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 Fig. 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 the FRW metric in the radiation and matter dominated eras using the PRS algorithm. The evolution a Minkowksi simulation is presented in Fig. 2 along with the corresponding spatial variation of R µ R µ , which would be zero in the case of a neutral vacuum solution, v + ≡ 0. In these simulations domain walls do not form in the individual fields of the linear representation. Rather, these fields evolve such that walls/condensates are manifest in the components of R µ . Snapshots of the spatial variation of each of the components of R µ during the same (2+1) dimensional simulation are presented in Fig. 3. Domains are mapped using the sign of R 1 as one would expect given the kink profiles of Fig. 1. The domain walls are manifest in the components of R µ due to their invariance under a local unitary transformation of the doublets. 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 Fig. 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 Fig. 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 [47,53], which has a common heavy Higgs mass,  [47,53], we assume here a Type-I Yukawa realization, so that the stringent constraints, mainly due to B-meson observables, can be avoided [54,55]. For such a scenario, the domains and spatial variation of R µ R µ are presented in Fig. 5.
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 Fig. 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 . Fig. 6 shows a a deviation from the t −1 scaling found in simpler models such as the Goldstone model [31]. 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 Figs. 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 [36]. 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.  For computational ease, the majority of our analysis has been performed in (2+1) dimensions. Nonetheless, we have also performed simulations of domain walls in (3+1) dimensions to verify that the same dynamical features are manifest in those simulations as well. A snapshot of a domain wall network from a (3+1) dimensional simulations of the Z 2 -symmetric 2HDM is given in Fig. 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 Fig. 8. We also find a deviation from the expected power law scaling as shown in Fig. 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 Figs. 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 and the corresponding force Using the equations of motion to substitute forΦ i andΦ † i in (4.23) and integrating, we find In order to evaluate the lower limit, x = b, we can derive asymptotic expressions for the kink and condensate. For a static kink/condensate solution in one dimension,Φ i (x) (with i = 1, 2), the Euler-Lagrange equations reduce to For a kink at x = a and anti-kink at x = −a we write the doublets as where δΦ 1,2 1 and Φ 0 1,2 are the CP-preserving vacua of (2.11). 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 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 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.32) Since the static solutions are asymptotically flat, their gradients vanish at infinity. Therefore, there is no contribution from the upper limit of (4.32). 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 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 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 Fig. 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 + . We also observe an increase in the energy per unit area as in Fig. 10 for solution 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 Fig. 12. In Figs. 11 and 12 we see that the vacuum values of the components R A remain unchanged, however, solutions with a relative transformation between the boundaries have features around the kink not observed in the minimum energy solution of Fig. 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 Fig. 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 Figs. 14 and 15, respectively. The profiles of the vacuum manifold parameters are consistent with those obtained for the kinks with   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 rotated vacua. Specifically, the peak in v + and the profile of ξ seen in Fig. 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 Fig. 16. We see that there is some clockwise winding of the electroweak angle, θ, around the domain wall. Note that Fig. 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 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, g H i ZZ = cos βO 1i + sin βO 2i , (5.12) 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.
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 [56] 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 Fig. 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. For both systems of equations a kink forms in R 2 and components of R µ have asymptotic values given by (2.23). The quantity R 2 is odd under a CP1 transformation whilst R 0 , R 1 and R 3 are even. Hence, we may define the gauge-invariant topological current, and corresponding topological charge,

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 Fig. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 18, the corresponding physical observables for the parameters chosen in Fig. 19  2D simulation of the evolution of domain walls in the CP1-symmetric 2HDM with Minkowski metric.
Parameters are {µ 2 1 , µ 2 2 , m 2 12 , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , λ 7 } = {1, 1, 0.1, 1, 1, 2, 1.5, 2, 0, 0}. Simulation was run for time, t = 480 with temporal grid spacing, ∆t = 0.2 and spatial grid size, P = 1024 with spacing, ∆x = 0.9. Plots progress in time left-toright and each plot is at double the timestep of the previous. The upper figures present the sign of R 2 and the lower ones show the spatial variation of R µ R µ . Clearly R µ R µ = 0 on the walls where R 2 changes sign. and FRW metrics in radiation and matter dominated eras are presented in Fig. 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 .

CP2 Symmetry
A CP2 transformation of the Higgs doublets is given by [30,45,46] Using the parameter constraints of Table 1, the CP2-symmetric 2HDM potential can be written as where we have introduced the notation Re (λ i ) = R i and Im (λ i ) = I i . The minimization conditions [16] for the CP2-symmetric potential are and 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 Re λ 6 e iξ v 2 1 + 3v 2 2 cot β , (6.7) 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 [49]. 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 Fig. 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 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 Fig. 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 Fig. 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 Fig. 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 Fig. 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 spon- taneous 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 the 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 FRW metrics, we have observed a deviation from the t −1 power law scaling found in simpler models. Specifically, in all three cases, 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 [36]. 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. 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.