Domain walls in the Two-Higgs-Doublet Model and their charge and CP-violating interactions with Standard Model fermions

Discrete symmetries play an important role in several extensions of the Standard Model (SM) of particle physics. For instance, in order to avoid flavor changing neutral currents, a discrete $Z_2$ symmetry is imposed on the Two-Higgs-Doublet Model (2HDM). This can lead to the formation of domain walls (DW) as the $Z_2$ symmetry gets spontaneously broken during electroweak symmetry breaking in the early universe and domain walls form between regions whose vacua are related by the discrete symmetry. Due to this simultaneous spontaneous breaking of both the discrete symmetry and the electroweak symmetry, the vacuum manifold consists of two disconnected 3-spheres. Such a non-trivial disconnected vacuum manifold leads to several choices for the vacua at two adjacent regions, in contrast to models where only the discrete symmetry gets spontaneously broken and the vacuum manifold consists of several disconnected points. Due to this, we end up with several classes of DW solutions having different properties localized inside the wall, such as charge and/or CP violating vacua. We discuss the properties of these different classes of DW solutions as well as the interaction of SM fermions with such topological defects leading to different exotic phenomena such as, for example, the top quark being transmitted or reflected off the wall as a bottom quark.


Introduction
Domain walls are a type of topological defects that arise after the spontaneous symmetry breaking (SSB) of a discrete symmetry. The formation of topological defects on a cosmological scale after SSB was already hypothesized in the last century by Zel'dovich [1], Kibble [2] and Zurek [3]. In the case that such a spontaneous symmetry breaking occurred in the early universe, causally disconnected regions may end up with different vacua located in the different sectors of the vacuum manifold of the theory. The boundaries between these regions are then called domain walls, which are sheet-like two-dimensionful surfaces where the discrete symmetry is restored. The vacuum manifold of the Standard Model (SM), for instance, is a 3-sphere which is topologically trivial. Therefore, no topological defects arise and beyond standard model physics with extra symmetries are needed in order to produce topological defects at the early universe.
Extensions of the SM scalar sector may lead to the existence of topological defects such as cosmic strings and domain walls. Such defects appear for example in axion models [4][5][6][7], supersymmetric models in the context of the spontaneous breaking of the R symmetry [8] and the Z 3 discrete symmetry in the NMSSM [9], discrete flavor symmetries [10,11] as well as in the Two-Higgs-Doublet model (2HDM) [12][13][14]. In this work we investigate the domain walls solutions in the 2HDM, where the SM scalar sector is extended with a second Higgs doublet of SU(2) L and charged under U Y (1). Such an extension allows several discrete symmetries in the model [12], whose spontaneous breaking after the electroweak symmetry breaking (EWSB) at the early universe can lead to the formation of domain walls. Recently, it was found in [13,14] that the domain walls in the 2HDM can have non trivial structures inside them. In particular, it was demonstrated that one dimensional domain wall solutions (usually denoted as kink solutions in the literature) exhibit CP and charge-violating vacua inside the defect. The spontaneous breaking of SU L (2) × U(1) Y alongside the discrete symmetry Z 2 leads to a degeneracy in the choice of the boundary conditions that one can impose on the vacua of different domains. This will then lead to several classes of kink solutions with different internal structures [15]. Such effects were already investigated for domain walls solutions arising in Grand Unified Theories such as SU(5) × Z 2 [15][16][17].
In this work we expand the analysis done in [14] to also include the variation of all Goldstone and hypercharge angles of the SU(2) L × U(1) Y symmetry, as well as to study the evolution of the 1D kink solution when using von Neumann boundary conditions. We then discuss the dependence of the kink solutions on the physical parameters of the 2HDM such as the masses of the extra Higgs bosons and the parameter tan(β). We also study the scattering of SM fermions off the different types of domain walls. This is done by solving the Dirac equation on a domain wall background that includes charge and CP-violating vacua, demonstrating the CP violating scattering of the fermions off the wall. We finally discuss some interesting phenomena such as top quarks being transmitted through the wall or reflected off the wall as a bottom quark. Such effects are charge violating and the charge difference will be absorbed by gauge fields living on the wall.
Our paper is organized as follows: First we briefly introduce the 2HDM, its potential and its different types of vacua. In section 3 we study the kink solutions in the model. We discuss in detail the different types of kink solutions and the effects of the variation of hypercharge and goldstone modes on the properties of the defect's core. We also discuss the profile of the kink solutions for different parameter points of the model. In section 4 we discuss the scattering of fermions off these kink solutions and show that the scattering exhibits CP and charge violating phenomena. Our summary and conclusions are given in section 5.

The Two-Higgs-Doublet Model
In this section we briefly introduce the general 2HDM and the used notation in this work. In the 2HDM the Standard Model Higgs sector is extended by an extra doublet charged under SU(2) L × U(1) Y . The general renormalizable scalar potential invariant under the SM symmetries is given by: Depending on the choice of the parameters, the potential can also be invariant under various discrete or continuous symmetries relating the Higgs doublets Φ 1 and Φ 2 [12]. The general Yukawa sector of the theory is then given by [18]: where ψ i and ψ j denote the different fermion generations. It is generally not possible to diagonalize the Yukawa couplings of fermions when having them couple to both Higgs doublets, this will then lead to the Yukawa couplings y 1 ij and y 2 ij being not simultaneously diagonalizable. Therefore, couplings between quarks of different flavor are then possible and this will lead to flavor changing neutral currents (FCNCs) at tree level [18]. Such phenomena are, however, strongly constrained experimentally, which makes it necessary to forbid them in the 2HDM. To avoid this problem one can impose that fermions with the same quantum numbers couple to the same Higgs doublet while other fermions couple to the second one. This can be achieved by imposing a Z 2 discrete symmetry on the Yukawa sector according to which the scalar doublets tranform in this way: There are therefore 4 types of 2HDMs depending on the choice of scalar doublets that fermions couple to [19]: Table 1: Types of Yukawa couplings between the fermions and the scalars in the 2HDM and the charges of the fermions under the Z 2 symmetry [18]. Q and L denote left-handed quark and lepton SU (2) L doublets while u R , d R and l R denote SU (2) L right-handed singlets.
The scalar potential which respects SU(2) L × U(1) Y × Z 2 is then given by: After electroweak symmetry breaking, the Higgs doublets acquire a vacuum expectation value. The 2HDM includes 8 scalar degrees of freedom. In our work we adopt the non-linear representation [13,14] to parameterize the vacua: where U is an element of the SU(2) L × U(1) Y group that is given by: where θ,g 1,2,3 are respectively the hypercharge angle and Goldstone modes, σ i denote the Pauli matrices, the generators of SU (2) and v sm is the vacuum expectation value of the SM Higgs doublet. Using this representation, we can separate the Goldstone and hypercharge modes from the physical vacuum parameters v 1 , v 2 , v + and ξ.
There are 3 possible types of vacua in the 2HDM: charge-breaking, CP-breaking and neutral: • The most general one occurs when v + is non zero and the vacuum is therefore charge breaking:Φ Such a vacuum is obviously non-physical as it will give the photon a mass and allow charge breaking interactions such as the decay electrons into neutrinos or the decay of top quarks into bottom quarks via new decay channels [20].
• The CP-breaking vacuum occurs when the phase ξ between the two doublets is non zero:Φ (2.9) this type of vacuum leads to CP-breaking Yukawa couplings and can be useful in the context of baryogenesis to generate the needed CP-violation.
In [21,22], it was shown that if a parameter point leads to a neutral minimum, then such a minimum of the potential will always lie above any possible charge or CP breaking minima. Throughout this work we will only consider that all regions of the universe ended up with a neutral vacuum after electroweak spontaneous symmetry breaking.  shows the dimensionless potentialV 2HDM = V 2HDM /(m 2 h v 2 sm ) with rescaled parametersv i = v i /v sm . One distinguishes degenerate vacua that can be related by a Z 2 transformation, as for example, (v 1 , −v 2 ) and (v 1 ,v 2 ) as well as degenerate minima that are related by a hypercharge transformation U Y (1) such as (−v 1 , −v 2 ) and (v 1 ,v 2 ). We also have a multitude of other degenerate vacua that can be obtained from the latter ones, by using a gauge transformation of SU(2) × U Y (1). The particle content of the CP-conserving 2HDM includes 5 physical Higgs scalars: two CP-even with masses m h and m H , one CP-odd with a mass m A and two charged Higgs bosons with a degenerate mass m C . Using the mass parametrization, one can trade the parameters in the scalar potential with physical parameters such as the masses of the physical scalars, the ratio between the 2 vevs of the doublets tan(β) = v 2 v 1 , the standard model vev v sm = 246 GeV and the mixing angle α. The potential parameters are therefore given by: (2.14) The term m 2 12 softly breaks the Z 2 symmetry leading the formed domain walls to be unstable and therefore annihilate some time after their formation. In this work we set m 2 12 = 0 and leave the effects of a small non-vanishing value for future studies.

Domain walls in the 2HDM
After electroweak symmetry breaking, the Higgs doublets acquire a neutral vev (2.10). This vacuum breaks SU(2) L × U(1) Y × Z 2 into U(1) em . In this case the vacuum manifold of the theory is homeomorphic to the coset space: which is topologically equivalent to two disconnected 3-spheres [12] as depicted in Figure  2a. The vacuum manifold has then two disconnected sectors related by a Z 2 symmetry. These sectors are non-trivial (in contrast to the case when only Z 2 is broken) and consist of vacua which are related by SU(2) L × U(1) Y transformations. This leads to the formation of different classes of domain walls due to the multiple choices that can be taken for the electroweak matrix U inside the two regions (see Figure 2a), in contrast to standard Z 2 domain wall solutions where the choice of the vacua inside the two domains is fixed to be vacua that are only related by a discrete symmetry [23]. Figure 2b shows different possibilities for the boundary conditions after EWSB. In order to get a kink solution to the scalar field configuration, the boundary conditions at ±∞ need to lie on disconnected sectors of the vacuum manifold. Starting at x → −∞ with a vacuum Φ − on the vacuum manifold sector M − corresponding to one 3-sphere: we end up at x → +∞ with a vacuum Φ + on the vacuum manifold M + . Fixing our choice for Φ − , we have multiple choices for the vacuum Φ + :  where U is an element of the broken electroweak symmetry group SU(2) L × U(1) Y and (v * 1 ) 2 + (v * 2 ) 2 = 246 GeV. In order to compute the kink solution of the field configuration interpolating between those two vacua, we need to minimize the energy of such a field configuration. In the 2HDM, the energy functional is given by: The first two terms describe the kinetic energy of the vacuum configuration, while the third is the potential of the scalar sector. There is an interplay between the kinetic energy that arises due to a changing profile of the field configuration as a function of x and the potential energy of this field configuration. Using the non-linear representation for the Higgs doublets (defined in (2.5), (2.6) and (2.7)) in the expression of the energy functional (3.4) we end up with: where, in terms of the vacuum manifold parameters: Writing down U (x) in terms of the Pauli matrices, we get: (3.7) where: In [14], the choice of the matrix U (x) was simplified to the case where only a single Goldstone mode or the hypercharge angle θ was allowed to be non-zero and have asymmetric boundary conditions at ±∞. We will expand the results to the general case, where all modes in U (x) can change. This will lead to more effects inside the domain walls compared to [14].
Since the calculation of all terms in (3.5) is straightforward but lengthy, we give the final expression of the energy functional in terms of the vacuum parametrization: where, 16) and the prime symbol denotes derivatives with respect to x. In order to get non-trivial vacuum configurations corresponding to different vacua at ±∞, we need to minimize the space integral of this energy functional with respect to small deviations in the fields: For static solutions, this leads to a system of differential equations analogous to the equations of motion: where ϕ i denotes the 8 fields in the two doublets. Since the energy functional and the equations of motion in the general case are very complicated, it is simpler to explain the behavior of the different fields inside the domain wall by choosing some special cases which make the energy functional and equations of motion considerably simpler. In the following sections we will consider a few special cases for the choice of U (x) at the boundaries. These simplified cases capture the interesting effects that influence the fields inside the domain wall. We will first start by considering the case of standard domain walls where the matrix U (x) = I 2 , constant everywhere in space, this means that we only consider vacua related by the Z 2 symmetry. The cases where only one mode of U (x) is different in the two domains while taking all the other modes to be 0 was already considered in [14]. Here we expand those results by first considering the case where more than one mode of U (x) changes across the domains and later the general case where the hypercharge and the Goldstone modes are chosen arbitrarily. We also discuss the behavior of the fields that interact with the kink solution for v 2 inside the wall.
To solve the 8 differential equations describing the vacuum configuration as a function of x, we use the same numerical algorithm used in [12][13][14], namely the gradient flow method. The gradient flow method introduces a fictitious time parameter to the vacuum profiles ϕ i (x, t). We then modify the minimization condition (3.17) to: Using this method, we find the solution of the system of differential equations by iteratively minimizing the energy functional at each time for the given boundary conditions until the vacuum configuration ϕ i (x, t) leads to a minimum in the energy. As we approach such a minimal energy configuration, the derivative of ϕ i (x, t) with respect to the fictitious time approaches zero and therefore we recover the static equations of motion in (3.17). The solution is then declared as found if, after several iterations, the vacuum configuration stays the same. We also adopt the same rescaling of the dimensionful parameters that was used in [12][13][14]: where m h = 125GeV denotes the mass of the SM Higgs particle. Such a rescaling is useful to get dimensionless space variablex and numerical values of the order 1.

Standard Domain Wall Solution
We start with the standard kink solution, where the two domains have the same and constant angles θ and g i while the vacuum expectation value of v 2 changes sign. The kink solution interpolates between the following two vacua: (3.21) We use (3.9) and set the derivatives for θ and g i to zero. We obtain the following energy functional that has to be minimized: (3.22) leading to the following equations of motion for the field profiles: As E(x) is independent of θ and g i , they therefore remain constant. The system of differential equations is solved in an interval −10 <x < 10. In Figure 3 we   to a positive value in the region on the right while crossing the value 0. This is the usual behavior of a kink solution. Note that also the value of the vacuum expectation value of v 1 (x) is also affected as both v 1 and v 2 are coupled via the differential equations.
To understand this change in v 1 (x) inside the domain wall of v 2 (see Figures 3a and 3b), we derive the standard energy functional E standard (x) in 3.22 taking the derivatives of all fields other than v 1 and v 2 to be 0 as they all vanish at all x: The effective mass term for v 1 is then given by: Inside the core of the wall (x = 0), v 2 = 0 and M 1 (0) = m 2 11 /2. When plotting the effective mass M 1 as function of x, we can see that the value becomes less negative inside the domain wall for both parameter points (see Figures 4a and 4b). Consequently the scalar potential V 2HDM (x = 0) = V 0 (v 1 , v + ) inside the domain wall will have its minima for v 1 at a smaller value and the potential barrier between the minima will be lower. Note that for other parameter points, (λ 3 + λ 4 + λ 5 ) can be positive and therefore the effective mass M 1 (0) inside the domain wall will be more negative than at the asymptotic values. This then corresponds to a bigger value for the minimum v 1 (0) of V 0 (v 1 ) and therefore v 1 inside the domain wall would be bigger than its asymptotic values at ±∞. However, a changing profile for v 1 (x) inside the wall will always lead to a positive contribution to the kinetic part of E(x). Therefore, one needs to make sure that this solution for v 1 (x) is stable. This is done by considering small fluctuationsṽ 1 (x, t) around the background kink    solution v 1 (x). The equation of motion describing such fluctuations is given by: Taking small fluctuations up to first order, this reduces to: for a fluctuation of the formṽ 1 (x, t) = e iwtṽ 1 (x). If a solution with w 2 < 0 exists, then the fluctuation v 1 (x, t) ∝ ew t ,where w = iw, grows with time leading to the instability of the solution for v 1 (x). For w 2 > 0 the fluctuation keeps oscillating around the found solution for v 1 (x). The behavior ofṽ 1 (x, t) for PP I is shown in Figure 5a, showing an oscillating fluctuation around 0, which means that the solution is stable. Note that a vacuum configuration, where v 1 (x) = v 1 is constant, will be unstable and the field v 1 (x) For the fields v + (x) and ξ(x), we observe that they stay equal to zero everywhere. A non-zero phase ξ provides a positive contribution to the energy functional (3.22) leading to a higher energy solution. In other words, ξ(x) = 0 presents the lowest energy solution.
Concerning v + (x), the situation is more complicated: in order to study its behavior inside the wall, we consider the terms in E that depend on v + (x), (cf. 3.22): The effective mass term for v + in the background of the DW is given by: Outside the domain wall, the effective mass is obviously positive and the potential minimum for v + is 0. Inside the wall, one gets which, depending on the parameter point can also be negative (see Figures 6a and 6b). The field v + can develop a condensate inside the wall, i.e. when M + (0) < 0. In that case the potential V 0 (v 1 (0), v + ) will have two non zero minima for v + (see Figure 6c). If such parameter points exist, it is energetically more favorable for v + to get a condensate inside the wall. However, the kinetic energy part is minimized for v + (x) = 0. Therefore, one should make sure that the kinetic contribution due to the spatial derivative of v + is   showsM + for PP II. Here the effective mass is positive everywhere and no charge violating vacua are expected inside the wall; (c) shows the potential inside the DW as a function of v + for PP I. In this case, the minimum is non-zero; (d) the same potential for PP II. In this case the global minimum of v + is zero.
not too large so that the solution for the non-zero condensate is stable inside the wall. In order to investigate the stability of the v + (x) = 0 solution inside the wall, we consider the linearized time-dependent equation of motion for a small fluctuationṽ + (x, t) around v + (x) in the background of the domain wall: For a small fluctuation of the formṽ the evolution of the fluctuation follows the differential equation: The small fluctuation around v + = 0 is unstable if w 2 < 0, making v + (x, t) ∝ ew t , where w = iw, growing with time and leading to the instability of the solution v + (x) = 0. In case w 2 > 0, the fluctuationṽ + (x, t) oscillates around the stable solution. Figure 6a shows one parameter point where M + (0) is negative inside the wall. However, the lowest energy solution has v + (x) = 0, meaning that the kinetic energy contribution from a v + (x) condensate inside the wall leads to a higher total energy than the solution with a vanishing v + (x) overall. We also verify numerically (see Figure 7) that the frequency w of a fluctuatioñ

Variation of a single angle across the wall
We now consider the effects of the Goldstone and hypercharge modes on the domain wall solution. We start by following the same approach in [14] and simplify U(x) (2.7) by allowing the variation across the wall of either the hypercharge angle θ(x) or a single Goldstone mode g i (x) at a time. In this case, the vacua at ±∞ will be rotated relative to each other, by either a U Y (1) transformation or a transformation related to one Goldstone mode of SU (2) L . For each case, we will discuss the solution of the equations of motion using either Dirichlet or von Neumann boundary conditions, where the former keeps the vacua at the boundaries fixed while the latter allows the dynamical variation of the vacua at the boundaries but keeps the spatial derivative of the vacua to be zero on the boundaries. For the following numerical solutions, we work in the alignment limit and fix the parameter point: m H = 800 GeV, m A = 500 GeV, m C = 400 GeV, tan(β) = 0.85 and α − β = 0. (3.37)

Variation of hypercharge θ
We first discuss the variation of the hypercharge angle θ(x) across the wall. In this case, the matrix U(x) (2.7) is given by: which leads to the energy functional E(x):  . We observe a non zero phase ξ(x) inside the wall when using Dirichlet boundary condition. This means that the vacua inside the wall are CP-violating (see (2.9)).
One can see immediately from (3.38), that a change in the hypercharge across the wall cannot lead to a charge breaking solution v + (0) ̸ = 0 inside the wall. The term 1 2 v 2 + (x)(dθ/dx) 2 in (3.38) always leads to a positive contribution to the energy and therefore it only minimizes the energy of the vacuum configuration when v + (0) = 0. Using the equation of motion for the hypercharge θ(x), one can derive a relation between the change in the hypercharge θ(x) and the derivative of the CP-violating phase ξ(x) [14]: Such an equation is only valid for finite energy solutions where the spatial derivatives of the vacua at the boundaries vanish. From (3.39) one would expect that, as the hypercharge angle starts to change from the value it has in one domain to the value in the other domain, the value of the phase ξ(x) will become non zero. Figure 8a shows the numerical solution of the equations of motion using Dirichlet boundary conditions with θ(−∞) = 0 and θ(+∞) = π/2. The initial guess for the profile of θ(x) was taken to be a hyperbolic tangent function interpolating between 0 and π/2. The numerical results point to a nonvanishing phase ξ(x) localized inside the wall, where the hypercharge changes significantly. However, this vacuum configuration has a higher energy than the vacuum configuration of the standard domain walls with ξ(x) = 0 (see Figure 9b). Therefore, for these parameter points, such domain walls are unstable and should decay into the standard domain wall configuration. One can also notice that the derivative in θ(x) outside the wall are non zero due to the choice of Dirichlet boundary conditions. This shows that minimizing the energy of the wall tends to change the hypercharge angle at the boundaries. For this reason, it is more advantageous to use von Neumann boundary conditions when solving the differential equations, so that the hypercharge angle in each domain can change dynamically to minimize the energy of the wall. This can be seen in Figure 8b, where the hypercharge in both regions evolve to be equal to each other. This will then lead to dθ/dx = 0 and therefore to a vanishing ξ(x) for all x. In Figure 9a we plot the profile of the solution at an intermediate time step using von Neumann boundary conditions. This shows how the hypercharge angle θ(x) at the boundaries change dynamically to minimize the energy and we see that the value of the CP-violating phase ξ(x) inside the wall gets smaller. Even though the CP-violating domain wall solution is unstable, it is expected that after EWSB, the hypercharge angles on causally disconnected regions of the universe can be different. Therefore the early stages of the formation of the domain wall network would exhibit such CP-violating vacua inside the wall until the profile of the hypercharge angle θ(x) relaxes to a solution where it is constant for all x.
We will also see later that in the realistic scenario, where we consider U(x) to be a general SU(2)×U Y (1) matrix, the stable domain wall solution will exhibit a non zero (albeit small) CP-violating vacua because the hypercharge angle will be different on both domains.

Variation of g 1
We now discuss the case when only the Goldstone mode g 1 is non-zero and different on both domains. The matrix U (x) (2.7) in such a case is given by: The corresponding energy functional E (3.9) is simplified to: In [14], by using the equation of motion for g 1 (x), an expression relating how the change in g 1 will affect v + and ξ was derived, namely: This would then imply that a change in g 1 across the wall will lead to a non-zero v + (x) and ξ(x) inside the wall.  The numerical solution of the equations of motion using Dirichlet boundary conditions, however, points to a vanishing v + (x) and ξ(x) (see Figure 10a), such a behavior was also found in [14]. This can be attributed to the fact that using Dirichlet boundary conditions for this system of differential equations is not the correct approach. The equation (3.42) is valid for a static field configuration that satisfies the equations of motion, has a vanishing derivative for the Goldstone mode at the boundaries and minimizes the energy functional dx E(x). However it is clear that the condition of vanishing derivative at the boundaries is not satisfied for our solution using Dirichlet boundary conditions, even though the solution is static (dg 1 /dt = 0). In Figure 11, we compare the energies of domain wall solutions using different boundary conditions. We observe that the energy of the solution using the Dirichlet boundary conditions is only a local minimum and has a higher energy than the standard domain wall solution. When using von Neumann boundary conditions (see Figure 10b), the Goldstone modes g 1 in both domains change dynamically to become the same value, which eventually leads to dg 1 /dx = 0. Such a field configuration is the correct solution: it explains the vanishing values of v + (x) and ξ(x) inside the wall, satisfies the relation (3.42) and minimizes the energy.  Figure 11: Evolution of the energy of the vacuum configuration with iteration time for Dirichlet and von Neumann boundary conditions. For g 1 having asymmetric boundaries, the energy is higher than the energy of the standard vacuum configuration.

Variation of g 2
We now discuss the case when only g 2 changes across the wall. The matrix U (x) (2.7) is given by: The energy functional E (3.9) simplifies to: One major difference in comparison with the analysis of the standard domain wall solution in section 3.1 is that dg 2 /dx ̸ = 0 induces a linear term for v + (x) in the potential (3.44): which can give a negative contribution to the energy of the vacuum configuration depending on the sign of dg 2 /dx and dv 2 /dx. Using the equation of motion for g 2 (x), one derives an expression relating the change in g 2 (x) to the derivative of v + inside the wall [14]: .
This implies that a variation in the Goldstone mode g 2 (x) will lead to a non-vanishing v + (x). In this case, even if ξ(x) = 0, it is possible to get a negative contribution to the energy of the wall by having a non vanishing v + (0), in contrast to the previous case where the dependence was on sin(ξ). A non-zero derivative for g 2 (x) leads to the creation of a stable v + (x) condensate inside the wall, if the energy E of such a solution is smaller than the energy of the standard domain wall solution.   Figure 12 shows the numerical solutions to the equations of motion using Dirichlet (Fig.12a) and von Neumann boundary conditions (Fig.12b). Like in the previous cases, we take the boundaries of the Goldstone mode to be: g 2 (−∞) = 0 and g 2 (+∞) = π/2. For the initial . Using Dirichlet boundary conditions, the relation is not fulfilled because the used boundaries at ±∞ give an unstable solution. Using von Neumann boundary condition, we get a perfect agreement.
guess, we use a hyperbolic tangent interpolating between both values. We observe a nonvanishing value for v + (x) inside the wall. The U (1) em is therefore broken inside the wall leading to exotic phenomena such as charge breaking processes and the photon getting a mass [24,25]. When using Dirichlet boundary conditions we observe again that the spatial derivative of g 2 (x) outside the domain wall is not vanishing. We can see in Figure 13a that the relation (3.45) is not exactly fulfilled using these boundary conditions, reflecting the instability of the solution with Dirichlet boundary conditions. Using von Neumann boundary conditions, we can verify that the non-vanishing v + (x) inside the wall is stable and that the change in g 2 (x) between the two regions gets enhanced, leading to a slightly higher value for v + (0) inside the wall compared with the solution using Dirichlet boundary conditions. We also verified that the relation (3.45) is satisfied for this choice of boundary condition (see Figure 13b). Even though the choice of von Neumann boundary condition is the correct choice to get stable solutions for our system of differential equations, we can still use Dirichlet boundary condition in order to get an approximate solution for a fixed choice of g 2 between the two regions, as this type of boundary condition fixes those values. In the case of von Neumann boundary conditions, if we choose different initial values for g 2 on the boundaries, we always end up with a solution for g 2 (x) where the difference ∆(g 2 ) stable = g 2 (+∞) − g 2 (−∞) is fixed. This means that when starting with random ∆(g 2 ) initial between [0, 2π], ∆(g 2 ) always relaxes to a fixed value ∆(g 2 ) stable that only depends on the mass parameters and tan(β). For the used parameter point (3.37), the v + (x) condensate is stable and the vacuum configuration has a lower energy than the standard domain wall solution (see Figure 14a). There are also other parameter points where the contributions from the derivative dg 2 /dx and a non-vanishing condensate v + (x) leads to a higher energy than the standard domain wall solution with v + (x) = 0. In such a   case, the Goldstone mode g 2 dynamically changes its values until it becomes equal in both domains leading to the charge breaking domain wall to decay into the standard domain wall solution as the latter has a lower energy (see Figure 14b).  Figure 15: Numerical solution of the DW equations of motion in the case of variation of g 2 . In contrast to the previous case, the derivative dg 2 /dx is negative, leading to a negative condensate v + (x) inside the wall.
For the case when g 2 (x) decreases when going from the region v 2 < 0 to the region v 2 > 0, for example, when taking g 2 (−∞) = π/2 and g 2 (+∞) = 0 (see Figure 15), we obtain a negative value for the condensate v + (x).

Variation of g 3
We now discuss the case when we allow g 3 to have different values on the boundaries. The matrix U (x) (2.7) is given by: The energy functional E (3.9) simplifies to: The change in g 3 (x) will lead to a change in the phase ξ(x) as was derived in [14] using the equation of motion for g 3 (x):  boundary conditions. We observe that the phase ξ(x) is non zero inside the wall, which means that the vacuum inside the domain wall is CP-violating. Using Dirichlet boundary conditions, we see again that the derivative of the Goldstone mode g 3 (x) is non zero at the : Evolution of the rescaled energy of the vacuum configuration using different boundary conditions. Notice that the energy using Dirichlet boundary conditions is higher than the energy using von Neumann boundary conditions. The class of domain walls that are CP-violating due to a variation in g 3 decays after some time to the class of standard domain walls.
boundaries: such a CP-breaking solution has a higher energy than the standard domain wall solution (see Figure 17) and is therefore unstable for this parameter point. The gradient flow method gives us a solution where the Goldstone field tries to change its value at the boundaries, reflecting the instability of the solution. However, if we use von Neumann boundary conditions, the Goldstone mode g 3 dynamically changes its value at the boundaries in such a way that, after some time, both regions end up having the same g 3 (see Figure 16b). The CP-breaking solution at the wall will then decay and we end up with a standard domain wall. Nevertheless, the Dirichlet boundary condition provides a good approximation for determining the amount of CP-violation inside the wall just after the formation of the defect.

Variation of the hypercharge angle θ and the Goldstone modes
We now turn to the case where we allow the hypercharge angle θ of U (1) Y to vary on both domains alongside the Goldstone modes of SU (2) L . In this case, we expect the formation of domain walls which are both CP and charge violating at the same time. We first discuss the case when the hypercharge angle θ and one single Goldstone mode g i are different on both domains. We then discuss the case of pure SU (2) L and finish with considering the general case when U (x) is a general electroweak symmetry matrix of SU (2) L × U (1) Y . We provide the numerical solutions for the equations of motion using both Dirichlet and von Neumann boundary conditions.

Variation of θ and g 1
We start by considering the case when the vacua of the two domains have different hypercharge angle θ and Goldstone mode g 1 . The matrix U (x) (2.7) simplifies to: Such a case is relevant to see whether a solution ξ(x) ̸ = 0 inside the wall will lead to a non-vanishing v + (x) condensate on the wall. The energy functional (3.9) is given by: (3.50) Figure 18a shows the numerical results for the vacuum configuration using Dirichlet boundary conditions. We take the values of the hypercharge angle θ and Goldstone mode g 1 to be: θ(−∞) = 0, g 1 (−∞) = 0 for the vacuum Φ − on the left and θ(+∞) = π/2, g 1 (+∞) = π/2 for the vacuum Φ + on the right. In this case we find a very small non-vanishing condensate v + (x) inside the wall. We also observe that the CP-violating phase ξ(x) is slightly enhanced inside the wall compared to the case where only the hypercharge varies across the wall (cf. Figure 8a). However, this CP and charge breaking solution is unstable as it has a higher energy than the standard domain wall solution.     Figure 19: Solutions for the equations of motion using von Neumann Boundary conditions. Notice that the hypercharge angle and the Goldstone mode g 1 are the same on both domains.
When using von Neumann boundary conditions (see Figure 19), both the hypercharge angle θ and the Goldstone mode g 1 evolve to become equal on both domains and the solution does not exhibit CP or charge violation.

Variation of θ and g 2
We now vary both θ and g 2 . The matrix U (x) (2.7) simplifies to: U (x) = e iθ(x) cos g 2 (x)/2 sin g 2 (x)/2 − sin g 2 (x)/2 cos g 2 (x)/2 . (3.51) In this case one would expect that we get both effects of charge and CP-violation inside the wall. The energy functional (3.9) is: (3.52) Figure 20 shows the solution for both Dirichlet (in Figure 20a) and von Neumann (Figure 20b) boundary conditions. The initial boundary conditions for θ and g 2 (θ(−∞) = g 2 (−∞) = 0, θ(+∞) = g 2 (+∞) = π/2) lead to a charge and CP-breaking vacuum inside the wall when using Dirichlet boundary conditions. However, such a solution is energetically unstable since a non-zero ξ(x) gives a positive contribution to the energy of the domain wall. We also observe a non-zero derivative for θ(x) and g 2 (x) at the boundaries.   This behavior reflects the instability of this solution, since the values of θ and g 2 on both domains try to change in order to further minimize the energy of the vacuum configuration. Using von Neumann boundary conditions (20b), where the boundaries can change dynamically to minimize the energy of the vacuum configuration, the derivatives of θ(x) and g 2 (x) at the boundaries vanish. One notices that for the solution which minimizes the energy the most, the charge breaking vacuum inside the wall gets enhanced while the CP-breaking phase ξ(x) inside the wall will start decreasing and eventually vanishes once the values for θ on both domains become equal to each other (as is shown in Figures 20d  and 20c). This vacuum configuration is stable, since its dimensionless energyÊ = 0.476 is lower than the standard domain wall's energyÊ standard = 0.507. This means that the standard domain wall solution will decay into the charge-breaking domain wall solution as can be seen in Figure 21. The Goldstone mode g 2 changes dynamically from a vacuum configuration where it is equal on both domains to a lower energy configuration with different values for g 2 on both domains. The hypercharge angle θ, however, stays zero on both domains and a CP-violating phase ξ(x) does not develop inside the wall, as such a solution would otherwise give a positive contribution to the energy of the defect.  For other parameter points, such as PP II (3.27), we observe the opposite behavior: the energy of the vacuum configuration for the CP and charge-breaking wall is higher than the energy of the standard wall. Such a scenario is shown in Figure 22 using Dirichlet ( Figure  22a) and von Neumann (Figure 22b) boundary conditions. The Dirichlet solution exhibits a CP and charge-breaking vacuum which is energetically unstable. Using von Neumann boundary conditions, the values of v + (0) and ξ(x) decrease until they vanish (both regions end up with the same value for g 2 and θ) and we recover the standard domain wall solution.

Variation of θ and g 3
We now discuss the variation of θ alongside g 3 . In this case one expects only a CP-violating solution and no charge violation inside the wall. The matrix U (x) (2.7) simplifies to: The energy functional (3.9) is given by: Using the equations of motion for θ(x) and g 3 (x), we can derive an equation that describes how the change in θ(x) and g 3 (x) causes a change in ξ(x): From such a relation, one expects the possibility of having an interference in the contributions of θ and g 3 . Figure 23 shows the numerical solution to the equations of motion. The initial condition were taken to be θ(−∞) = g 3 (−∞) = 0, θ(+∞) = g 3 (+∞) = π/2. In this case, we again see that the Dirichlet boundary condition leads to a localized CP-violating phase ξ(x) inside the wall and that the derivatives of the Goldstone modes do not vanish at the boundaries. This vacuum field configuration is energetically unstable and decays to the standard vacuum configuration with ξ(x) = 0 for all x. In contrast to the results from the variation of θ and g 3 individually, the minimum vacuum configuration gets two different values for θ and g 3 in the two domains (see Figure 23b).   Using (3.55), one sees that if the derivative of θ(x) is equal to half the derivative of g 3 (x), one obtains dξ/dx = 0. This clarifies why the values of ξ(x) are vanishingly small despite the different Goldstone modes on both domains (see Figure 23b). This condition can be verified in Figure 24, where we see that both expressions agree numerically.

Variation of SU
(2) L modes g 1 , g 2 and g 3 In this special case the vacua in both regions are related by an SU(2) L gauge transformation and we ignore the effects coming from the change in the hypercharge angle θ. The matrix U (x) (2.7) simplifies to: (3.56) The energy functional (3.9) is given by:  The numerical results using von Neumann boundary conditions are shown in Figure 25. The initial profile of the Goldstone modes were chosen to be g i (−∞) = 0 and g i (+∞) = π/2, (where i denotes 1,2,3) with a tangent hyperbolic function interpolating both boundaries. The solution has a stable charge-violating vacuum v + (x) inside the domain wall, as well as a small and stable CP-violating phase ξ(x). In contrast to the previous results, ξ(x) in this particular case is asymmetric. This behavior could be attributed to the fact that the profiles of the Goldstone modes g i (x) are not symmetrical (see Figure 26). There is also an interference between the Goldstone modes as some have negative derivatives while the others have a positive derivative inside the wall. Note that, in this case, the Goldstone modes keep having a different value on both domains.  Figure 26: Goldstone modes g 1 (x), g 2 (x) and g 3 (x). Notice that the profiles of the Goldstone modes are asymmetric inside the wall. This suggests an interference between the different Goldstone modes.

General Domain Wall Solution
We finish the discussion of domain wall solutions in the 2HDM by considering the case with a general matrix U (x) (2.7): We recall the general formula for the energy functional 3.9: Here, the Goldstone modes g i and hypercharge angle θ are chosen randomly and can be different on both domains.
Boundary at +∞ positive positive 0 0 π/2 0 π/6 0 Boundary at −∞ positive negative 0 0 0 0 π/6 0 Table 2: Asymptotic values of the fields at the boundaries. The initial profile for θ(x) is taken to be a tangent hyperbolic function interpolating between 0 at −∞ and π/2 at +∞. We use von Neumann boundary condition to get the lowest energy solution. Table 2 shows an example of the initial asymptotic values for the fields at ±∞ that we use for the numerical calculations with θ having a tangent hyperbolic profile interpolating the values on the two boundaries. We use von Neumann boundary conditions to get the numerical solution of the 8 equations of motions describing the profiles of the fields.

Figures 27 and 28
show the numerical solution of the vacuum configuration for this choice of hypercharge angle θ and Goldstone modes. The solution features a stable charge-violation as well as a small but stable CP-violating phase ξ(x) inside the wall. In contrast to all previous cases, the behavior of the hypercharge angle θ(x) and the Goldstone modes g 1 (x) and g 3 (x) is non-trivial inside the wall. We also note that, even though we started with the Goldstone modes g i being the same on both domains, the lowest energy solution has different values for g i (±∞).      One notices the non-trivial behavior of the hypercharge angle and Goldstone modes g 1 (x) and g 3 (x).

Dependence of the kink solution on the parameter points of the 2HDM
In this subsection, we briefly describe the dependence of the kink solution on the masses of the Higgs bosons of the 2HDM as well as tan(β). For all the discussed parameters, we use the alignment limit: α = β. In this work we do not take into account experimental constraints on the 2HDM. Our choice of the parameter points discussed is done in order to give an overview of the different properties that arise for domain walls in the 2HDM.
We start by analyzing the effects of the parameter points on the standard domain wall solution. We therefore take the matrix U (x) to be the identity.  We observe that for higher m H , the value of v 1 inside the domain wall (v 1 (x = 0)) becomes smaller. This is explained by the behavior of the effective mass M 1 (x) inside the domain wall (see 3.29) as shown in Figure 29c. Consider a potential of the form: Outside the wall and for parameter points with masses m H > 80 GeV, the effective mass term M 1 is smaller (more negative) than inside the domain wall. Therefore, the minimum v DW

1
= v 1 (x = 0) of the potential V 0 (v 1 , v 2 = 0) inside the wall (depicted in blue) is smaller than v * 1 , the minimum of the potential V (v 1 , v 2 = v * 2 ) outside the wall (depicted in orange) as is shown in Figure 30a. Note that for m H = 80 GeV, the opposite effect occurs and v 1 (x = 0) inside the wall is bigger than outside of it (see Figure 30b). This is due to the fact that for m H = 80 GeV, M 1 (x) becomes more negative inside the wall than outside of it (see Figure 29c). This leads the minimum v DW

1
= v 1 (x = 0) of the potential V 0 (v 1 , v 2 = 0) inside the wall to be bigger than v * 1 , the minimum of the potential V (v 1 , v 2 = v * 2 ) outside the wall as is shown in Figure 30b. Notice also that the effective mass term M 1 (x = 0) is the same for all mass parameters. This is because the parameter m 11 does not depend on m H . As for the dependence of the kink solution for v 2 (x) on m H , we observe that increasing the mass leads to a thinner profile for the kink as can be seen in Figure 29b. Figure 29d shows M + (x) for different m H . A negative effective mass M + (0) inside the wall and therefore the possibility of having a charged condensate v + (x) localized on the wall is possible for higher masses of H.
We now vary the mass of the charged Higgs m C while keeping all other masses fixed to 200 GeV and tan(β) = 0.85. We do not observe any change in the profiles of the vacua v 1 (x) and v 2 (x) (see Figure 31a and 31b). However, we observe a change in the effective mass M + (x) for different masses m C and it becomes smaller inside the wall as shown in Figure 31c. However, the valueM + (x = 0) is only negative for small masses of the charged Higgs. Therefore, if the mass m H of the CP-even Higgs is small, it is possible to have a stable charged condensate in the wall only if the charged Higgs masses m C is very low.   We also study the effect of m A and m H on the CP-violating phase ξ(x). In this case we study the scenario when the two domains have different hypercharge angle θ and use the Dirichlet boundary condition with θ(−∞) = 0 and θ(+∞) = π/2. We first fix the masses m H = 800 GeV, m C = 400 GeV and tan(β) = 0.85 and vary m A . The results are shown in Figure 33a and we observe that increasing m A leads to a smaller phase ξ(x) inside the DW. However, when fixing m A = 500 GeV and varying m H , we observe the opposite behavior: increasing m H leads to a higher CP-violating phase ξ(x) (see Figure 33b).      Figure 34a). This leads the scalar potential inside the wall V 0 (v 1 , v + ) to have its minimum at v + (0) = 0 (see Figure 6d). This means that any charge violating solution for those domain walls is unstable. In order to get a stable v + (x) condensate inside the wall for this low value of tan(β), we need to choose high values for m H . As we increase tan(β), the fraction of parameter points with a negativê M + (0) increases (as is shown in Figure 34b) and a charged condensate inside the domain wall can be stable (in case of different g 2 Goldstone mode on both domains). For high values of tan(β), most parameter points have a negativeM + (0) as can be seen in Figures  34d and 34c. In the case that v + develops a stable condensate inside the wall, the masses m H and m C can have a sizable effect on the maximum value v + (x = 0) of such a condensate. In order to study this, we solved the equations of motion for the case when the vacua have different values of g 2 using von Neumann boundary conditions. This was done for two scenarios: • To study the effect of m H , we fix m C = 400 GeV, tan(β) = 0.85 and vary m H between 80 GeV and 1100 GeV. The results are shown in Figure 35a. We observe that the v + condensate is unstable for m H < 580 GeV. The value of v + (0) inside the wall increases with the mass m H .
• To study the effect of of m C , we fix m H = 800 GeV, tan(β) = 0.85 and vary m C between 80 GeV and 1100 GeV. The results are shown in Figure 35b. In this case, v + (0) gets smaller with heavier m C . For m C > 680GeV the condensate becomes unstable.

Non-Topological kink solutions in the 2HDM
We discussed in the previous sections the possibility of having topological kink solutions that interpolate between vacua belonging to disconnected sectors of the vacuum manifold. We now briefly discuss the scenario when both vacua at x → ±∞ belong to the same sector (in the case of the 2HDM, the same 3-sphere M + or M − ) as shown in Figure 36. We consider the following two degenerate vacua: (3.62) One can obtain Φ − by performing a U (1) Y transformation on Φ + . Therefore, both these vacua belong to the same sector M + . A kink solution for such a vacuum configuration is not topologically protected against a variation in the fields and therefore any such solution should, in principle, be unstable [26]. Indeed we obtain a kink solution (see Figure 37a) for this configuration using von Neumann boundary conditions. The solution has a higher energy than the topological kink solution, where the vacua are related by a Z 2 transformation (see Figure 37b). However the solution seems to be stable, even when using von Neumann boundary conditions. This might be explained by the fact that a potential energy barrier  exists between the vacua Φ + and Φ − (see Figure 1a). If we transform the vacuum Φ + with an SU(2) × U(1) Y gauge transformation U (2.7), the kink solution can develop charge and/or CP-violating vacua localized inside the wall but such a vacuum configuration is no longer stable. In the 1D simulation, one vacuum grows in the region of the other vacuum. The charge and/or CP-violating vacua also move alongside the wall. Therefore, after some time, the kink solution dynamically decays into the trivial vacuum configuration, where one vacuum occupies the whole space.

SM fermions scattering off the domain walls
The existence of different types of domain walls in the 2HDM can have profound implication on the physics of the early universe. For instance, it was shown in [25] that photons with small frequencies such as in the CMB will scatter off superconducting domain walls with a large reflection coefficient. In [23], it was shown that domain walls in Grand Unified Theories arising after the spontaneous breaking of SU(5)×Z 2 , interact with the Higgs scalar field. This interaction induces exotic scattering phenomena of fermions off the domain wall via the Yukawa sector, such as neutrinos being reflected as d-quarks, with the electric and color charges being absorbed by gauge fields living on the wall. In our work, we study the interactions of SM fermions with 2HDM domain walls via the Yukawa sector by considering the scattering solution of the Dirac equation of SM fermions within a domain wall background. We leave the case of fermion zero modes and bound states solutions [15] on the walls, which might be relevant for baryogenesis scenarios, for future work. The solution for the scattering of fermions off standard domain walls generated by the spontaneous breaking of a discrete symmetry can be found in [23,27]. One finds that for thin walls, the rate of reflection and transmission of fermions off the walls is: In this work, we want to get analytical solutions for the Dirac equations within the background of different types of domain walls arising in the 2HDM. As the functions describing the spatial dependence of these vacuum configurations are non-trivial, it is appropriate to use a thin-wall approximation to simplify the form of the different vacuum configurations inside the wall. The thin-wall approximation is valid for the scattering of fermions which have wavelengths larger than the width of the wall. As another simplification, we do not consider the back-reaction effects of fermions on the vacuum configurations in our study, which could change the spatial kink profile of the vacuum configuration [28].

Thin-Wall approximation
In this section, we briefly discuss the validity of this approximation for domain walls in the 2HDM as well as describe the vacuum configuration in this approximation. From results of the last chapter we can infer a typical width of the domain walls to be approximately L w ≈ 4/m h . The typical wavelength λ f of a particle in the thermal plasma in the early universe is proportional to 2π/T , where T is the temperature of the thermal plasma. The thin wall approximation is valid when L w < λ f corresponding to particle momenta smaller than 200 GeV, which is sufficient to describe the momenta of most particles existing in the thermal plasma after EWSB. Taking a thin-wall profile for the domain walls, we can approximate the kink solution of v 2 (x) to be a step-function: As for the vacua v 1 (x), v + (x) and Im(v 2 (x)e iξ(x) ), it is possible to approximate them with a delta distribution: v 2 (x) describes the profile of the imaginary part of the CP-violating mass term that appears in the Dirac equations (see Figure 38). We saw in the last chapter that the phase ξ(x) vanishes (in the simplified cases) when using von Neumann boundary conditions (see e.g sections [3.2.1] and [3.2.4]). However, the vanishing of ξ(x) inside the wall only occurs after the dynamical evolution of the Goldstone modes or hypercharge angle θ inside the different domains, making them equal to each other. Therefore, the CP-violating phase should still be substantial at the moment of the formation of the domain wall network and the study of CP-violating scattering of fermions off the wall is relevant for that period, when considering the simplified cases. However, as shown in sections [3.3.4] and [3.4], for domains rotated by a relative SU (2) L or an electroweak symmetry SU (2) L × U Y (1) the CP-violating phase ξ(x) of the kink solution is small but stable. For simplicity, we consider the simplified cases in this work and leave the scattering of fermions off general domain walls in the 2HDM for future investigations.

CP-Violating interactions of fermions with the domain walls
We start with the case of fermions scattering off CP-violating domain walls induced by a difference in the hypercharge angle θ. Recall that in this case, the Higgs doublets are given by: where the matrix U (x) is given by: We can remove this matrix U (x) from the Yukawa sector by performing a gauge transformation, leading to a pure gauge term for the hypercharge gauge field B µ : The Yukawa Lagrangian for the up-type quarks in the type-2 2HDM is then given by: where Y u,L and Y u,R are the hypercharges of the left and right-handed up-type quarks respectively, y u is the Yukawa coupling of the Higgs doublet to the up-type quark and u L , u R are the left-handed and right-handed components of the up-type quark respectively. One can then derive the Dirac equation for up-type quarks: where P L and P R are the left and right-handed projector operators and m R (x), m I (x) are the real and imaginary parts of the fermion mass, respectively: 14) where we apply a small angle approximation for the CP-violating angle ξ(x), which is motivated by the simulations of the domain walls in the previous chapter. After multiplying the left hand side of (4.12) with iγ 1 , we can rearrange this equation and get: The solution to this equation can be calculated in an analogous way to the case of fermion scattering off a CP-violating bubble wall [29]: whereP is an ordering operator. We consider the case of a plane wave solution of a quark scattering off the domain wall from the left (x < 0), that can either be reflected or transmitted to the other region (x > 0): u(t, x) = e −iEt+ipux u tra for x > 0, (4.19) where E denotes the energy of the incoming quark and u inc , u ref and u tra are 4-component spinors describing respectively the incident, reflected and transmitted fermion. (4.20) By plugging this ansatz into the Dirac equation (4.12) for the regions far away from the wall (where the vacua v 1 (x) and v 2 (x) take on their asymptotic values and ξ(x) = θ(x) = 0), we can derive relations between the different components of the spinors: (4.21) The matching condition of the solution u(x) at x = 0 is calculated using (4.17): where ϵ is a small number taken to ϵ → 0. The final result is given by: where: (4.24) Note that the spinor u(x) is not continuous at x = 0. This is due to the delta-functions iñ v 2 (x) (see (4.5) and Figure 38). It is a known issue that the presence of delta distributions in the Dirac equation leads to a discontinuity of the spinor's wave function [30][31][32][33] (in contrast to the discontinuity in the derivative of the wave function when dealing with deltadistribution potentials in the Schrödinger equation). In the limit of a standard domain wall, a → 0 andÂ is a zero matrix, therefore we recover the continuity condition of the quark's spinor at x = 0 given by u(−ϵ) = u(+ϵ). The Dirac spinor of an incident particle moving in the positive x-direction is given by: (4.25) Using (4.25) alongside the equations (4.21) and (4.23), we can find the solution for the spinor components: , with: In order to get the transmission and reflection coefficients of particles scattering off the wall, we need to calculate the fermion currents on both regions: Using the expression (4.25) for the incident spinor, we can derive the transmission and reflection coefficients for the up-type quark scattering off the wall: In Figure 39, we show a comparison between the reflectionR and transmissionT coefficients of quarks scattering off the CP-violating DW and the coefficients R and T for quarks scattering off standard CP-conserving DW for two different cases. The first case depicted in Figure 39a shows the results for a 1 = 2.1 and b 2 = 1, while for the second case (shown in Figure 39b) the parameters are: a 1 = 2.1 and b 2 = 3. We see that in both cases, the reflection and transmission coefficients differ a lot from the standard reflection and transmission coefficients. In particular, we see that the reflection coefficient for incident particles with higher momenta grows and stays the dominant process. Since the analytical results for the general case are quite complicated, we consider some special cases in order to understand the physical interpretation of the solution. We first consider the case where b 2 is small compared to a 1 . In such a scenario, the reflection and transmission coefficients simplify to:R As shown in Figure 40a, in this case and for high values of a 1 , the reflection coefficient increases with momentum. This seems counter-intuitive as particles reacting with a potential barrier should have a higher probability of crossing when they have higher energies. One can also deduce from (4.36) that for big values of CP-violation (a 1 is big), the reflection coefficient approaches 1 for all momenta, as the terms proportional to cosh 2 (2a 1 ) will be dominant.
(b)T as a function of a 1 and momenta Figure 40: Results for the reflection (a) and transmission (b) coefficients for the case when the numerical value of the Yukawa contribution a 1 is much higher than the numerical value from the gauge field contribution b 2 . Notice that for higher values of CP-violation a 1 , the reflection coefficient increases with momentum whileT decreases.
For the opposite case b 2 >> 2ṽ 2 , the transmission and reflection coefficients are: , (4.38) .  For this case, the reflection and transmission coefficients will oscillate with a b 2 dependence (see Figure 41). Notice that for b 2 = π all particles will be transmitted, irrespective of their incoming momentum. In Figure 42, we show the oscillatory behavior of the reflection and transmission coefficients for this special case as a function of b 2 for different momenta of the incident particle. We see that for small momenta (see Figure 42a), the deviations of the scattering rates from those of the standard DW are only relevant for b 2 in the vicinity of b 2 = π. In contrast with particles with higher momenta (see Figure 42b).
We now calculate the rate of reflection and transmission of the left handed (LH) and right handed (RH) components, derived from the currents of left and right-handed particles: leading to the transmission and reflection coefficients: , (4.47) where: From equations (4.26)-(4.29) it is clear that for a 1 = 0 (as in the case of the standard domain wall), (u * 1p u 2p ) and (u * 1m u 2m ) are purely imaginary expressions leading to the vanishing of the second terms in the equations for T L,R and R L,R . In such a case the reflection and transmission coefficient do not depend on the chirality of the particle. Because analytical expressions for the general case are lengthy and complicated, we give for simplicity the analytical expressions in the limit when b 2 → 0 : Note that for small momenta, these rates can become negative (see Figure 43a for a 1 = 0.5). This is not an issue as these rates are not describing probability amplitudes for the particles and onlyR andT should be, in principle, positive numbers between 0 and 1. In case a 1 → 0, we do not observe CP-violation in the scattering rates of the particle. In order to get a rate for the CP-violation, we calculate the difference ∆ LR = R L − R R . The results are shown in Figure 43b. We find from (4.49)-(4.52), that when the momentum of the incoming particles gets larger, the left and right-handed rates should converge to the same value. This explains the behavior of ∆ LR becoming smaller at higher momenta. Since we get different rates for R L and R R , the motion of the wall will generate an axial asymmetry in front of it.

Charge violating interactions of fermions with the domain walls
We now discuss the case when the domain wall exhibits a non negligible value for v + inside its core. In this case the vacuum inside the domain wall is charge-breaking and the Yukawa sector inside the wall now includes couplings between fermions with different charges. In the following, we consider the case of a Yukawa sector in the type-2 2HDM. The relevant terms are: where Φ c 2 = iσ 2 Φ 2 . y u , y d and y l denote the Yukawa couplings for the up-type quarks, down-type quarks and leptons respectively. Q L and L L denote the left-handed quark and lepton doublets charged under SU L (2) and finally u R , d R and e R denote the right-handed up-type, down-type quark and lepton singlet under SU L (2). We take the simplification that only one of the Goldstone modes changes between the two domains. In this case, it is the second generator of SU(2) which leads to charge violating effects inside the domain wall and the CP-violating phase ξ(x) vanishes (see 12b). The Higgs doublets are: The matrix U (x) in this case is given by: We can remove this matrix from the Yukawa sector by performing a gauge transformation, leading to a pure gauge term for the gauge field W µ 2 confined inside the wall: Writing the Lagrangian in terms of the individual quark fields, we get: We then derive the Dirac equation for the up and down-type quarks: To solve this system of equations, we first re-write it in the matrix form introduced in the previous chapter, cf. (4.16): where:Ĝ Notice that all terms that mix the up and down-type quark only include the left-handed component of the down-type quark. This will induce a chiral asymmetry as we will discuss later. We consider the scattering of an incident up-type quark off the domain wall. In the standard case, when v + = 0, the particle can be either reflected or transmitted as an up-type quark. However, due to the mixing between the up and down-type quarks in the Dirac equations (4.60) and (4.61), there is now the possibility that the incoming up-type quark gets also reflected or transmitted as a down-type quark after its interaction with the domain wall, inducing a charge violating interaction and the difference in charge between the up and down-type quark will be absorbed by the gauge bosons living inside the wall [15,23]. Using a plane wave solution for the spinor fields, and inserting these expressions in (4.62), we get: (4.69) which can be solved by taking: whereP is, as in the previous chapter, taken to be an ordering operator. We now derive the expressions for the spinors in the different regions x < 0 and x > 0. Far from the wall, the Dirac equations are: In the region x < 0, we parameterize (4.71) and (4.72) by: leading to the following relations between the components of the spinors: where p u denotes the momentum of the incident up-type quark and p d denotes the momentum of the produced down-type quark. For the region x > 0, we get: i / ∂u + y u v 2 u = 0, (4.78) leading to: We obtain for the complete spinor: where we used a plane wave solution for the incident up-type quark moving in the positive x-direction: Since the analytical formulas in the general case are very complicated and lengthy, we present the results in Figure 44a for the numerical values of m = 0.5, k = 10 and a mass m u = 172.76 GeV for the top quark and m d = 4.2 GeV for the bottom quark. We observe that the reflection and transmission probabilities as a bottom quark are non-zero. In Figure  44b, we can see that as we increase the momentum of the incoming particle, the rate of top quarks being transformed into bottom quarks after the interaction with the domain wall becomes higher, while the probability of the quark staying a top quark decreases.   We also verify that all the reflected or transmitted bottom quarks are left-handed as it can be already deduced from equations (4.60) and (4.61), where the coupling to the bottom quark includes only the left handed projector on the spinor. As can be seen in Figure 45 the rate of reflection and transmission coefficient as right-handed bottom quarks is zero. Concerning the chirality of the top quark after the interaction with the wall, we observe a difference in the transmission rate between the left-handed and right-handed components of the top quark, while the reflected top quarks do not show a chiral asymmetry (see Figures 45c and 45d). Therefore the scattering of top quarks on the charge breaking wall will generate a chiral asymmetry in front as well as behind the wall.   We observe that the reflection rate as a top quark does not depend on the chirality of the particle, while the transmission for the top quark is chirality-dependent.
For the case k >> m and m → 0 corresponding to having the numerical value of v + (0) inside the wall much bigger than the change in g 2 between the two domains, the formulas become considerably simpler. Using (4.25) for the spinor components of the incident top quark we get: where: The results for the transmission T t,b and reflection R t,b coefficients as top or bottom quarks  Notice that the reflection rate as a bottom quark is in this case higher than the transmission rate even for high momenta of the incoming incident top quark. We also observe that the rates are almost the same for k > 5.
are shown in Figures 46a and 46b. An interesting feature is that the rate of reflection as a bottom quark is higher than the rate of transmission as a bottom quark, even for higher momenta. Recall that we observed the same behavior in the previous section [4.2] when the numerical value of the imaginary massṽ 2 is dominant and the reflection rate for particles becomes higher than the transmission rate.
We now look at the case when the charge-breaking term k = y u v + is small compared with the change in the Goldstone modes between the two domains (m > k), cf.(4.89). The boundary condition at x = 0 is given by: One can see that, in this case, the transmission coefficient as a bottom quark gets higher as we increase the momentum of the incoming particle. The difference between this case and the one where the Yukawa term is dominant (k >> m) is that there is no chirality flipping of the particle in the coupling between the gauge fields living inside the wall and the fermions (see (4.59)). For m >> k one can get rather simple analytical expressions: T u (p u ) = −4 cos(m)(−1 + g 2 d ) 2 + (1 + 6g 2 d + g 4 d )(3 + cos(2m)) 8 cos 4 ( m 2 )[g d (d u + g u ) + (1 + g 2 d ) tan 2 ( m 2 )] 2 , (4.103) (4.104) In this case the rate of reflection and transmission coefficients will oscillate with increasing m. Such a behavior is shown in Figures 47c and 47d. In order to study this oscillating behavior in more detail, we fix the momentum of the incoming top quark and vary m between [0, 2π]. The results are shown in Figure 48. We observe that the rate of conversion of top quarks into bottom quarks also vanishes for m = π. We now consider the scattering of the second generation quarks off the wall. We take, as an example, the scattering of charm quarks off the domain wall. In this case the charge breaking parameter k gets smaller due to the small Yukawa coupling of the charm quark to the Higgs doublet. Figure 49 shows the results for k > m. In this case we observe that the charge breaking effect is very small and most charm quarks scatter off the wall as charm quarks, even for high values of v + = 65 GeV corresponding to k = 0.43. However, when we consider the case m > k, the charge breaking effect can be quite high depending on the value of m as is shown in Figure 50. Finally, we also mention that anti up-type quarks scattering off the domain walls will lead to exactly the same rate of transformation into anti down-type as the rate of the up-type quarks transformed into down-type quarks. Therefore, the interaction of fermions with these types of domain walls does not lead to a net generation of charge in the early universe, a strongly constrained phenomenon [34]. However, as we saw in the previous chapter, domain walls solution in 2HDM can exhibit simultaneously both a charge and a CP violation inside the wall. This will lead to the generation of a non-zero CP-violating phase ξ inside the domain wall. This phase, along with a non-zero v + inside the wall, might lead to the generation of a net charge in the early universe, as particles and antiparticles will interact with the wall with different rates. In such a case, the domain wall network has to annihilate very quickly in order to avoid generating a charge asymmetry higher than the observed cosmological constraints. Another possible problem with these types of domain walls is that they could efficiently deplete the number of up-type quarks into down-type quarks. For example, bottom and anti-bottom quarks could be generated from the interaction of top and anti-top quarks with the wall. Due to the large difference of the masses between the two flavors, the inverse reaction would be suppressed and we would end up with a surplus of bottom and antibottom quarks well before the usual annihilation time of top quarks in the thermal plasma. This phenomenon might then have consequences on Big Bang Nucleosynthesis (BBN).

Summary and conclusions
In this paper we investigated the different classes of domain walls arising in the 2HDM after EWSB. We extended the work done in [14] and included the variation of all the angles of the SU(2) L × U(1) Y symmetry. In contrast to the standard domain wall solution, where only a discrete symmetry such as Z 2 gets spontaneously broken, we saw that the breaking of abelian and non-abelian symmetries alongside the discrete symmetry leads to the formation of kink solutions with non-trivial effects in the core of the defect, such as CP and charge-violating vacua. We have found that these different classes of kink solutions can be unstable and decay to the standard kink solution if their energy is higher than the energy of the standard kink solution. We demonstrated this behavior using von Neumann boundary conditions, where the vacua, the hypercharge angle θ and the Goldstone modes g i can change their boundary values dynamically in order to minimize the energy of the field configuration. When discussing the simplified cases (where only the hypercharge or/and a single Goldstone mode g i were allowed to change on both domains) we found that for the CP-breaking kink solution, the CP-violating phase ξ(x) inside the defect decays after some time and we recover again the standard kink solutions. Nevertheless, such CP-violating effects can be quite sizable at the time of the domain walls formation as demonstrated when using Dirichlet boundary conditions instead of von Neumann conditions. In the case of charge breaking kink solutions, we showed that the stability of the solution depends on the sign of the effective mass M + of v + inside the wall. When investigating the general case where all the equations of motions of the Goldstone modes as well as hypercharge angle are taken into account, we saw in particular, that these modes had a non-trivial profile inside the wall in contrast to a kink-like profile that was found in all other simplified cases. We also found that the CP-violating phase ξ(x) in that case is stable even though it had a small value. We also focused on the interaction of Standard Model fermions with different classes of kink solutions. In the case of CP-violating kink solutions, the fermions scatter in a CP-violating manner on the wall. In particular, we observed that the transmission and reflection coefficients can deviate substantially from the ones found for the case of standard kink solutions. This is in contrast with results found for CP-violating bubble walls [35], where the transmission and reflection coefficients T and R do not depend on the CP-violating phase. However, in a realistic scenario, the phase ξ(x) is dynamically vanishing and the scattering of fermions off these CP-violating simplified kink solutions is only relevant at the first stage of the evolution of the domain wall network after its formation. The CP-violating phase eventually decays with time as the Goldstone modes and hypercharge angle θ change dynamically, becoming equal in both domains and the scattering of fermions off the walls eventually occurs in a standard way. Even though we showed that the CP-violating phase ξ(x) inside the wall (for the simplified cases) is unstable, the initial value just after the formation of the domain wall network can be large and the CP-violating scattering of fermions could have an important impact for the initial phases of the domain wall's evolution. This effect generates a chiral asymmetry in front of the wall, which could serve to achieve electroweak baryogenesis via domain walls. In fact, due to v 2 (x = 0) = 0 and the possibility of v 1 becoming smaller inside the wall, the sphaleron rate inside the wall is less suppressed than outside of it and the scenario of the annihilation of the domain wall network would satisfy the last Sakharov condition, i.e. the departure from the thermal equilibrium [36]. Electroweak baryogenesis using topological defects such as cosmic strings was already investigated in [37,38], where it was shown that the volume suppression factor due to the 1D nature of cosmic strings is a huge constraint. As domain walls are 2D sheets in space, the volume suppression factor is, in principle, much less constraining and studying electroweak baryogenesis using 2HDM domain walls is an interesting scenario. A model-independent investigation of such a possibility was done in [39,40] where the topological defects were formed before EWSB. Studying this formalism in the context of 2HDM, where the domain walls form just after EWSB is subject of future work. Furthermore, we considered the scattering of fermions off charge-violating kink solutions. We took, as an example, the scattering case of incident top quarks and showed that they can also be reflected or transmitted off the wall as bottom quarks. Such a scattering is charge-violating and the charge is absorbed by the gauge fields and bosons living on the wall. The rate of such charge-violation phenomenon is the same for the antiparticles. This means that we do not expect a global charge asymmetry from such phenomena. This might be different when considering the scattering off both CP and charge-breaking DW. As a charge asymmetry in the universe is heavily constrained [34], this might provide strong constraints on the 2HDM. Such a scattering would also, in principal, deplete the population of up-type quarks in the early universe, which might have some consequences on BBN. The consequences of such phenomena in the context of early universe cosmology are subjects of future work.