Domain walls in the extensions of the Standard Model
Abstract
Our main interest is the evolution of domain walls of the Higgs field in the early Universe. The aim of this paper is to understand how dynamics of Higgs domain walls could be influenced by yet unknown interactions from beyond the Standard Model. We assume that the Standard Model is valid up to certain, high, energy scale and use the framework of the effective field theory to describe physics below that scale. Performing numerical simulations with different values of the scale we are able to extend our previous analysis [1] and determine its range of validity.
We study domain walls interpolating between the physical electroweak vacuum and the vacuum appearing at very high field strengths. These domain walls could be formed from nonhomogeneous configurations of the Higgs field produced by quantum fluctuations during inflation or thermal fluctuations during reheating.
Our numerical simulations show that evolution of Higgs domain walls is rather insensitive to interactions beyond the Standard Model as long as masses of new particles are grater than . For lower values of the RG improved effective potential is strongly modified at field strengths crucial to the evolution of domain walls. For instance its minima become degenerate for around . We find that even in the case when the minima of the potential are nearly degenerate Higgs domain walls decayed shortly after their formation for generic initial conditions. On the other hand, in simulations with specifically chosen initial conditions Higgs domain walls can live longer and enter the scaling regime.
We also determine the energy spectrum of gravitational waves produced by decaying domain walls of the Higgs field. For generic initial field configurations the amplitude of the signal is too small to be observed in present and planned detectors.
a]Tomasz Krajewski, a]Zygmunt Lalak, b]Marek Lewicki, a]Paweł Olszewski
Prepared for submission to JCAP
Domain walls in the extensions of the Standard Model

Institute of Theoretical Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, Warsaw, Poland 
ARC Centre of Excellence for Particle Physics at the Terascale (CoEPP) & CSSM,
Department of Physics, University of Adelaide, South Australia 5005, Australia
Keywords: domain walls, gravitational waves, Higgs effective potential, nonrenormalizable operators, effective field theory
1 Introduction
The precise measurement of the Higgs boson mass at the Large Hadron Collider [2, 3] enabled study of the scalar effective potential in the Standard Model up to superplanckian field values. The potential at field strengths higher than approximately GeV drops below its value at the electroweak vacuum and develops a new minimum at superplanckian field strengths. This opens the possibility of tunneling from the physical electroweak vacuum to the deeper minimum at very large Higgs field value. This effect could render the electroweak vacuum unstable and has been widely discussed in the literature [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Metastability of the electroweak vacuum for present central values of measured parameters of the SM has been shown under the assumption that the SM is valid up to the Planck scale. However, the problem of stability of the electroweak vacuum can be modified, especially if particles with masses of the order of GeV or lower, not described by the SM, exist in nature. This scenario has been investigated in [6] in the framework of the effective field theory. Furthermore, even if the electroweak vacuum is metastable, existence of the second deeper minimum poses severe problems for cosmological models [16, 17, 18].
One of the most promising solutions to problems of the standard cosmology is the occurrence of an inflationary epoch. The analysis based on the FokkerPlanck equation [11, 19, 20, 17, 21] generally predicts that the Higgs field distribution after inflation is, in the first approximation, Gaussian with the standard deviation of the order of
(1.1) 
where is the value of the Hubble constant during the inflation and is the number of efolds which for experimentally viable models exceeds value of 50–60.
We will investigate the case in which Hubble constant during inflation was of the order of the position of the local maximum separating the two minima in the Higgs potential. Thus, the generated field fluctuations would be large enough to create a network of domain walls which interpolate between areas of the Universe occupied by the Higgs field lying in different minima of the potential.
We perform an analysis of the evolution of these structures using numerical lattice simulations based on the Press, Ryden, Spergel (PRS) algorithm [22].
We will refer to the basin of attraction of the electroweak vacuum, that is the range of fields whose absolute values are less than , as . Analogously will denote the range of fields of an absolute value greater than .
(1.2) 
Let us also adopt labels for the exact field values that minimize the potential: for the electroweak vacuum and for the high energy vacuum.
Domain walls, interpolating between nondegenerate minima are known to be generically unstable [23]. As we have shown in our previous article [1], assuming validity of the SM up to the Planck scale, networks of Higgs domain walls decayed in conformal time of the order of . Initial configurations of the Higgs field with at least 49 % of the expectation values belonging to decayed to the superplanckian minimum. Models of the early Universe predicting such configurations of the Higgs field are excluded by experiments. On the other hand for networks initialised with the higher fraction of lattice sites belonging to , the final state of their evolution was the electroweak vacuum. Moreover, the short decay time of these networks guarantee consistency with present experimental data.
In this article we investigate the possibility that there exist new, yet unknown physical phenomena which are not described by the SM. We parametrize the influence of these new interactions on the physics of the Higgs boson in the framework of the effective field theory (EFT). To that end we introduce a nonrenormalizable operator suppressed by the scale in the Lagrangian density of the SM. Performing simulations with different values of the initial standard deviation of the field strength and different values of the suppression scale , we are able to determine the influence of this new interaction on evolution of Higgs domain walls. We find that a naive expectation that new physics will not modify dynamics of Higgs domain walls if the suppression scale will be much higher then the position of the local maximum of the effective potential is accurate. Moreover we have considered lower suppression scales leading to a reduced difference between values of the effective potential in the two minima. The scenario with nearly degenerate minima is interesting, because it allows metastable networks of domain walls which would lead to a reach variety of phenomenological implications.
During the process of decay of domain walls, the energy of the field is transferred to other degrees of freedom and both SM particles and gravitational waves (GWs) can be produced. The next point in considering low suppression scales producing nearly degenerate minima is the possibility to increase the energy density of produced GWs. In our previous work we estimated the energy spectrum of GWs produced by decaying domain walls of the Higgs field in the SM valid up to the Planck scale. The obtained present energy density of GWs was orders of magnitude smaller than the sensitivity of present and planned detectors of GWs and so we excluded the possibility of their detection in the near future. However, our previous results do not apply to the case of metastable domain walls associated with nearly degenerate minima of the potential. Authors of [24] claim, basing on a semianalytical approximation [25], that GWs emitted from decaying Higgs domain walls in the case of nearly degenerate minima can be strong enough to be measured in planned detectors. In our numerical simulations we investigate evolution of Higgs domain walls with the potential with nearly degenerate minima. Unfortunately, we were not able to study the scenario proposed in [24] with very long decay time of the network of domain walls using lattice simulations due to their limited dynamical range.
Our simulations show that Higgs domain walls decay in conformal time shorter than after their formation if generic initial configurations of the Higgs field are assumed. We were however, able to simulate networks with decay times up to five times longer by setting very specific initial conditions (to compensate for the asymmetry of the potential around the local maximum) as well as a finetuned coupling constant for the interaction to obtain nearly degenerate minima. These examples show the possibility of formation of metastable networks of domain walls when the minima of the potential are nearly degenerate. Moreover, we prove that these longlived networks of domain walls evolve in the scaling regime, satisfying assumptions of the semianalytical method used in [24]. However, the formation of networks evolving in the scaling regime requires specifically chosen initial conditions in addition to assumptions used in [24]. Unfortunately, removing unphysical, short wavelengths modes in the initial configuration of the field what is necessary for computing GWs spectrum in finite lattice simulations, we also spoil the specific choice of the initial configuration required for the networks longevity.
The paper is organised as follows. In section 2 we briefly present the form of the RG improved effective potential of the SM in the presence of the nonrenormalizable operator . The influence of the modified potential on properties of Higgs domain walls is analysed in section 3. We discuss initial conditions for our simulations in section 4. Next, in section 5 we discuss dependence of the decay time of networks of domain walls on the suppression scale . Subsection 6 is devoted to an investigation of the possible metastability of networks of domain walls. Finally, we calculate the energy spectrum of GWs for the case of nearly degenerated minima in section 7. Our results are summarized in section 8.
2 RG improved SM effective potential with operator
In our previous studies [1] we have investigated evolution of domain walls of the SM Higgs field. The aim of this study is to understand the way in which our previous results can be modified by new, not yet discovered interactions of the Higgs boson.
The starting point of our considerations is the zero temperature scalar effective potential in the SM which takes the form [8]:
(2.1)  
(2.2) 
which depends implicitly on the running couplings , , , , , (only the largest Yukawa coupling of the third family is included). Explicit loop corrections to the quartic term are gathered in , where the field variable appears only as . Hence, the substitution allows one to resum large logarithms by solving renormalization group equations (RGEs) to appropriately high scales. We use threeloop RGEs and boundary conditions together with determined up to twoloops. This leaves us with the effective potential as a function of a single variable,
(2.3) 
This potential is famously unstable and exhibits a double nearcriticality in its dependence on the masses of the Higgs boson and the top quark, discussed for example in [8].
We consider modifications of that might stem from hypothetical existence of heavier particles, not included in the determination of the potential in the SM alone.
The EFT provides a computationally convenient framework for treating corrections from high energy physics to observables at lower energies via integrating out heavy states. This paradigm states that quantum corrections to scattering amplitudes of light modes at low energies coming from heavy states (whose masses are much higher relative to fourmomenta of the light modes) can be reconstructed (in a perturbative expansion) by inclusion of nonrenormalizable operators with properly adjusted coupling constants into the Lagrangian density of the light fields. Due to dimensional arguments, the coupling constants of the nonrenormalizable operators have to be suppressed by inverse powers of certain energy scale which is naturally of the order of masses of heavy states. Reversing the reasoning, the EFT approach can be used to parametrise, in a model independent way, the yet unknown high energy physics when treating values of coupling constants of nonrenormalizable operators as free parameters.
Basis of nonrenormalizable operators up to dimension 6 for SM fields was determined [26, 27] in the past. We will concentrate our analysis on the operator which influences the effective scalar potential the most since it is the only dimension 6 interaction term that simply enters the potential at the treelevel. Thus we choose to examine potential of the form:
(2.4) 
where is the energy scale used in the numerical integration of the RGEs [6]. We also define the scale which unambiguously parametrises the coupling constant of the nonrenormalizable operator . We simply use function as a oneparameter modification of with the robust characteristic that a significant deviation is switched on at a particular scale, . It is highly debatable, we realize, whether value of the coupling hypothetically determined in low energy expansion of nonrenormalizable corrections, valid for , can be meaningfully used in (2.4) to encode new physics in its entirety when is not much smaller than . Following other authors [28, 29, 30], we will still use as a free parameter which lets us smoothly interpolate between the scalar potential of the pure SM and the Higgs potential with nearly degenerate minima.
Scenario in which new interactions lift the value of the effective potential at the high field strength minimum up to the value at the electroweak one is supposed to have very interesting phenomenology. The general prediction states that networks of cosmological domain walls corresponding to a potential with degenerate minima should be metastable. Due to their long decay times, and consequently large abundance in late cosmological epochs, such networks might imprint sizable distortions in cosmological observables.
In the case of Higgs domain walls the most promising observable is the energy density spectrum of GWs produced by their decay. Previously we argued [1] that if the SM is valid up to the Planck scale, GWs produced by Higgs domain walls are orders of magnitude to weak to be detected in present or planned experiments. The main reason of such a small amplitude of the signal is the short lifetime of Higgs domain walls. Generally, the energy density of a network of domain walls decreases slower then the energy density of both: radiation and dust. Thus, the network’s contribution to the total energy density of the Universe grows in cosmic time. As a result, the longer the network lives, the higher abundance of GWs it produces.
The energy spectrum of GWs emitted by Higgs domain walls with nearly degenerate minima of the scalar potential has been previously studied in [24] with the semianalytical approach. This method is reliable if domain walls evolve in the so called scaling regime, during which the number of walls in the Hubble horizon is approximately preserved in time. The authors of [24] claim that GWs emitted from Higgs domain walls can be detected by future detectors if the difference between values of the Higgs potential in two minima is very small.
However, a small difference between values of the Higgs potential at minima requires a severe finetuning of the value of the coupling constant . The condition for nearly degenerate minima can be written down as:
(2.5) 
where is a small number that parametrizes the degeneracy of both minima and is a constant whose value is determined by the renormalization of the cosmological constant. After dividing both sides by and suppressing the constant we get:
(2.6) 
At this point the first conclusion from numerical studies of the RG improved potential of the SM with the operator is that nearly degenerate minima are obtained for suppression scales of the order of , so . Moreover the experimental data determines and . Thus .
The numerical study of the RG improved SM effective potential reveals that for nearly degenerate minima . This leads to the conclusion that a natural scale for is , so in order to produce nearly degenerate minima the coupling constant needs to be finetuned in order to get the value of the factor down to .
The difference between values of the RG improved potential at two minima as a function of the suppression scale is presented in the figure 1. The finetuning problem is made evident in this plot by the steep character of the presented function around the scale .
Finetuning of the value of coupling constant poses severe problems for numerical methods. Firstly, in order to calculate the value of the RG improved potential with finetuned numerical value of running coupling constants one is required to numerically solve RGEs of the SM with very high precision. Furthermore, the evolution of longlived networks of domain walls in the case of nearly degenerate minima is difficult to simulate numerically, because it requires a large number of integration steps leading to accumulation of truncation errors.
3 Properties of SM domain walls in the presence of the operator
We are interested in the evolution of the time and space dependent expectation value of the Higgs field . By definition, an expectation value of a field strength satisfies the equation of motion: , where is the 1PI effective action. The 1PI effective action of the SM is not known exactly and we approximated it by the following expression:
(3.1) 
where is a real scalar field which models the Higgs field in our simulations and is the effective potential of the SM in zero temperature with the nonrenormalizable operator included. Assuming the gravitational background in the form of the FriedmanRobertsonWalker metric background:
(3.2) 
where Latin indices correspond to spatial coordinates, is cosmic time and denotes conformal time (such that ), the eom for the approximated effective action from eq. (3.1) can be written down in the form:
(3.3) 
We have solved eq. (3.3) in numerical simulations on the 3D lattice, based on the PRS algorithm [22]. To properly model dynamics of a network of domain walls in these simulations, one needs to know the width of a typical domain wall . In the past, simulations with the physical width of walls varying from 2 to 100 lattice spacing (i.e. the physical distance between neighbouring points) were performed [22, 31, 32, 33, 23, 34, 35, 25].
In our previous article [1] we have proposed the algorithm for computing the width of domain walls for a generic potential. We started with the observation that the time independent, planar (i.e. translationally invariant in two directions) solution of eq. (3.3) in Minkowski background () is given in the implicit, integral form:
(3.4) 
Furthermore the potential energy density of the solution (3.4) is given by the following integral:
(3.5) 
Finally, we have found values and such that and the majority of the potential energy density of the solution (3.4) is stored between and :
(3.6) 
Our estimation of the width of domain walls is then given by:^{1}^{1}1We are using the convention , so lengths and times are expressed in units of .
(3.7) 
The estimated width of domain walls as a function of the suppression scale of the nonrenormalizable operator is presented in the figure 2. Resulting values of the width lay in the range .
We have chosen the physical lattice spacing to be equal to which leads to widths of walls contained in the range . Henceforth, we will often use as a unit of energy and inverse distance in spacetime.
4 Initial conditions for simulations
Inflationary models predict that the Higgs field expectation value during inflation can be statistically described by the probability distribution whose evolution satisfies the FokkerPlanck equation. Solutions of the FokkerPlanck for the Higgs field [31, 11, 19, 20, 17], in the first approximation, are of the form of the Gauss distribution:
(4.1) 
with the standard deviation of the order of
(4.2) 
where is the value of the Hubble constant during inflation and is the number of efolds. The Hubble constant is the main parameter of inflationary models. Solving majority of problems found in the standard cosmology requires at least from to efolds, so the standard deviation can be determined in an inflationary model. Unfortunately, evolution of turns out to be independent of the mean value of the Higgs field and this value had not changed significantly during inflation. As a result the mean value of the Higgs field after inflation cannot be determined in an inflationary model and a priori it can take any value from to the Planck scale.
We focus on the evolution of networks initialised with . This choice leads to the least conservative bounds on inflationary models. The lower the value of is, the higher the fraction of lattice points in is and consequently the more likely it is that the evolution of the network will end up in the electroweak vacuum.
Our simulations were initialised at the conformal time ranging from to . In our previous work [1] we showed that taking initial times below this interval does not modify the results. On the other hand, cosmological domain walls whose evolution we model in our numerical simulations, need to be superhorizon at the initialisation (i.e. width of domain walls must be larger than the Hubble horizon at the time of initialisation ). The conformal time when domain walls are formed in the early Universe, must be greater than for the the initial numerical fluctuations to be smoothed out by the evolution of the field.
5 Decay of domain walls
In our previous work [1] we investigated the evolution of networks of Higgs domain walls in the SM valid up to the Planck scale. We showed that obtaining a network that collapses to the electroweak vacuum is not as difficult as one may expect, although it requires certain tuning of initial conditions. Our simulations predict that only a slight dominance of the electroweak vacuum in the initial configuration is required in order for the evolution to end in this vacuum. Previously we concluded that only initial probability distributions (4.1) with lower then () produce networks that decay to the electroweak vacuum. Unfortunately, this value of is not a gaugeindependent quantity, however, the fraction of space occupied by the Higgs field strength belonging to is a gaugeindependent quantity [20]. Moreover, we found that networks of Higgs domain walls are highly unstable. The estimated time of their decay ranges from to . Such short lifetimes exclude a scenario in which SM domain walls dominate the Universe leading to large distortion of the Cosmic Microwave Background Radiation.
The main aim of this work is to generalize our previous results to include BSM physics. As we discussed in section 2, new interactions beyond the SM may significantly change the form of the Higgs effective potential. These modification influence the evolution of networks of Higgs domain walls. In this section we present results of our new numerical simulations. We study dependence of the decay time of networks on the suppression scale of the nonrenormalizable operator . An example of the time evolution of such a network is presented in the figure 3. In this figure the ratio of number of lattice sites occupied by field strengths in to the total volume of the lattice, as a function of conformal time , is plotted for two different values of the conformal initialisation time . Results calculated for various values of the suppression scale and the same initialisation conditions have been arranged. Presented results for each case were averaged over five simulations on the lattice of the size .
A few observations can be deduced from presented plots. Firstly, the evolution of domain walls quickly becomes insensitive to corrections from nonrenormalizable operators for the suppression scale larger than for which the minima of the potential are degenerate. This conclusion is supported by results of our scan over the space of the initialisation parameter .
Estimated dependence of the conformal decay time of the network on the value of the standard deviation at the initialisation and the value of the suppression scale is plotted in the figure 4. Blue regions in these plots were extrapolated from the simulations in which evolution of networks ended in the basin of the attraction of the EWSB vacuum and red ones from networks decaying into the basin of the attraction of the high field strength minimum (excluded by experiments).
Nearly vertical shape of resulting contours proves that dependence on the value of suppression scale of the decay time of networks is very weak in the shown range of . Moreover, a quick inspection of plots in the figure 4 reveals that late domain walls, i.e. initialised at , have the longer decay time what is consistent with our previous findings [1].
Secondly, from plots in the figure 3 we see that the decay of networks driven by the effective potential with nearly degenerate minima proceeds differently than it would with higher suppression scales . The initial value of the fraction in the case of is higher than for higher values of the suppression scale . Moreover, networks of domain walls in this case decay in the conformal time at least three times shorter. In the figure 5 we present lifetime of networks of domain walls driven by the effective potential with the nonrenormalizable operator suppressed by the scale close to the value at which the minima are degenerate.
First observation from the figure 5 is such that for the final state of evolution of networks initialised with vanishing mean vale is always the electroweak vacuum. Moreover decay time of networks decreases with the decreasing suppression scale . This behaviour disagrees with the prediction that domain walls should be metastable in the case of a potential with degenerate minima. This effect is the result of a change in the position of the local maximum which increases as decreases. Figure 6 illustrates how changes as a function of the suppression scale .
Second result of the shift in that can be observed in the figure 3 is a change in the starting value of the fraction in the case of . Figure 7 shows the decay time of networks of domain walls as a function of the suppression scale and the ratio . Again, the nearly vertical arrangement of contours indicates weak dependence on the suppression scale .
6 Metastable networks decay
Our simulations show that for generic initialisation probability distributions ( and different values of ) networks of Higgs domain walls were unstable and decayed shortly after the formation. Moreover, considering two values of the suppression scale finetuned in order to produce potentials with nearly degenerate minima with difference between values at minima of the order of and , we have not observed significant difference in the lifetime these networks.
Previous studies predict that networks of cosmological domain walls whose dynamics is driven by potentials with nearly degenerate minima evolved in the so called scaling regime. This means the number of domain walls in the Hubble horizon was preserved during the evolution of the network and the typical size of walls was of the order of the Hubble radius. These networks then became unstable due to a too large Hubble radius and their evolution ended by decay of the network. This behaviour was confirmed in many numerical lattice simulations [22, 36, 33, 37, 38, 23, 34, 35, 39, 40].
However, majority of simulations preformed in the past have dealt with spontaneous breaking of approximate global symmetries. In these scenarios the difference of values of potentials between their minima are generated by small corrections to symmetric potentials. Thus nearly degenerate minima implicate nearly symmetric potential in these situations. Obviously this is not the case of the SM with nonrenormalizable operators. The RG improved effective potential of the SM with the nonrenormalizable operator is asymmetric, even when minima are degenerate. We estimate the typical scale of the modulus of the derivative of the potential on both sides of the local maximum separating the minima (outside a close neighborhood of the maximum) to be of the order of:
where we neglected terms proportional to which are subleading for large field strength. The term in the bracket can be expressed via the RGEs in terms of running coupling constants of the SM, thus it depends logarithmically on far enough from its zero (i.e. the position of the maximum of ). The potential derivative behaves as , so that it grows with increasing field strength . As a result the RG improved effective potential is steeper on high field strengths side of the local maximum. This asymmetry of the potential is the source of the indicated discrepancy between our simulations of the evolution of Higgs domains walls and previously performed simulations of domain walls for spontaneously broken approximated global symmetries.
Previous numerical studies of dynamics of cosmological domain walls [23] revealed that asymmetry of the potential can be compensated by asymmetry of the initial configuration of the field leading to metastable networks. In order to validate this reasoning we have searched for asymmetric initial configurations producing longlived networks of Higgs domain walls for different values of the suppression scale . We have used probability distributions with and satisfying relation , where is the position of the local maximum, and is a free parameter. This parametrization of initialisation parameters guarantees that the fraction of lattice sites belonging to depends only on and is independent of (which depends on ). We performed scans with various values of and we found a significant increase in the decay time of networks for values higher than i.e. highly finetuned initial configurations.
Decay times obtained in our simulations initialised with are presented in the figure 8. Points with filled markers correspond to networks whose evolution ended in the electroweak vacuum and points with empty markers to ones for which the final state was the high field strength minimum. Points marked with blue circles corresponds to , with red squares to , green triangles are for and gray diamonds for . Each point in the figure 8 is an average over five independent simulations.
We observed large fluctuations of values of the decay time for longlived networks. This behaviour can be easily understood. The initial probability distribution is reproduced by the initial configuration on a finite lattice with finite accuracy. Small fluctuations of the initial configuration introduced by finite accuracy of the initialisation algorithm result in large relative fluctuations of the decay time in these finetuned scenarios.
Observed decay times were as long as and strongly point towards possibility of producing metastable networks for highly finetuned initial configurations. However, obtained long decay times need to be treated with caution, because the late evolution of networks in these simulations extends beyond the time range for which finite lattice simulations are reliable. For conformal times larger than the size of the lattice which for presented results was equal to , the modelled patch of the space corresponds to only a fraction of the Hubble horizon at the end of simulations. Thus, boundaries of the lattice are no longer causally disconnected and their influence on the late evolution of networks may be significant. Furthermore, even if boundary conditions would not disturb significantly the simulated evolution, results of long time simulations can still be inaccurate. In the scaling regime Hubble horizon contains only a few domain walls and it could happen that the simulated patch accommodates none of them. As a result the computed decay times may underestimate real lifetimes of longlived networks of Higgs domain walls.
For the value of equal to all considered networks decayed into the EWSB vacuum. For higher values of , i.e. lower fraction of field values belonging to at the lattice at the initialisation and more finetuned initial configurations, networks for certain higher values of suppression scales decayed to the (global) high field strength minimum. For majority of the considered values of the suppression scale the higher was the value of (the more finetuned was the initial configuration) the more stable were formed networks. Moreover, networks whose dynamics was driven by the effective potential which global minimum is EWSB minimum were less stable than networks for which the high field strength minimum is the global minimum of the effective potential.
However, for and the most stable networks were ones with the effective potential with nearly degenerate minima and the highly finetuned difference of the values at the minima. Moreover they were the most stable configurations observed in all our numerical simulations. The time evolution of these networks is presented and compared with the less finetuned potential in the figure 9. Dotted lines in the figure 9 correspond to simulations with the effective potential with the difference between values in minima and dashed lines to ones with values of the potential finedtuned to . Colours of lines indicate the value of in the probability distribution used to initialize simulations: blue for , red for , green for and gray for . Plot (a)a shows time dependence of the ratio of lattice sites occupied by the field on the electroweak side of the barrier to the volume of the lattice. Networks of domain walls with the potential with initialised with and decayed shortly after the formation. Networks initialised with higher values of , equal to and , decayed later with lifetime increasing for larger . Lifetimes of networks of domain walls with the potential with are longer than the ones with the less finedtuned potential initialised with the same values of . Networks initialised with were not very sensitive to the level of the finetuning of . However, evolution of networks initialised with parameters satisfying relation with differ strongly for two considered levels of finetuning. Networks initialised with driven by the potential with were the most stable networks observed in our numerical simulations.
Better insight into the evolution of longlived networks is provided by the plot (b)b which presents time dependence of the surface area of networks. From this plot one may deduce that the evolution of three fastest decaying groups of networks differs qualitatively from the evolution of the rest. The surface area of these three networks decreases monotically with the conformal time . On the other hand, in time dependence of surface area of the rest of simulated networks one finds a local maximum followed by the period of slow decay. Such behaviour was previously observed in lattice simulations [23] and was connected with networks entering the scaling regime. This behaviour of surface areas is the evidence for Higgs domain walls also evolving in the scaling regime.
A visualization of an exemplary network of Higgs domain walls and its evolution is shown in the figure 10. The isosurface of the field strength corresponding to the position of the maximum approximated by the interpolation from values at lattice sites is visualized for four different conformal times: and , , . In the first panel (a)a we see several domain walls spreading over whole lattice at the beginning of the scaling regime. The further decay of the network during the simulation, presented in panels (b)b and (c)c, ends with only one bubble of one vacuum in the background of the second vacuum (panel (d)d). The visualization was prepared with data from numerical simulations performed on a lattice with size . Percolation theory predicts that networks evolving in the scaling regime consist of walls spreading over the whole Universe. Thus, large sheets visible in the figure 10 shows that a metastable network can be formed for the specific initial configuration of the field.
Finally, we want to stress that even though this example of the evolution of Higgs domain walls in the scaling regime gives evidence for the occurrence of the generally predicted metastability of networks of domain walls with the potential with nearly degenerate minima, the realisation of this scenario in the early Universe was unlikely. As our numerical simulations show, the formation of the long lived network of Higgs domain walls required not only the Higgs effective potential with nearly degenerate minima, but also the specific initial configuration of the Higgs field. The problem of the production of the metastable network of Higgs domain walls corresponds to simultaneously finetuning of two parameters to unnaturally small values.
7 Gravitational waves
We use a numerical algorithm proposed in [41] for computation of the GWs produced during inflation and later used also in simulations of domain walls [34, 25, 1].
The algorithm is based on Einstein equation linearised around FRW metric background solution. We assume that domain walls evolve in the radiation domination epoch during which the scale factor grows linearly:
(7.1)  
(7.2) 
where is the value of the Hubble constant at the conformal time and is a constant which sets the value of the scale factor at the time when radiation starts to dominate the Universe. We neglect the constant in our numerical simulations, assuming that the evolution of domain walls took place deep in the radiation domination era, i.e. the ratio was small at the time when domain walls were formed.
We neglect the back reaction from domain walls to te evolution of the homogeneous gravitational background. This approximation is supported by the fact that the critical energy density in our simulations is much larger than the difference between values of the maximum and the electroweak minimum of the potential . We also assume that Higgs field fluctuations were the main source for GWs (constituted the dominant contribution to the transversetraceless part of the energymomentum tensor ). The lowest order in perturbations of metric around FRW background contribution to is given by
(7.3) 
where are projections on the transversetraceless part and is the field which models the Higgs field in our simulations.
The linearized Einstein equation can be solved using Green’s function and Fourier transform in the background given by eq. (7.2) if one neglects the back reaction on the Higgs field from GWs. The solution used in [41] is of the form
(7.4) 
where is the Fourier transform of with the normalization:
(7.5) 
and is the conformal time before which the source appeared. Moreover, in the Fourier space projection operators can be expressed in a more compact form:
(7.6) 
where .
The spectrum of gravitational waves’ energy density per unit logarithmic frequency interval can be calculated using [41]:
(7.7) 
with the function , computed in our simulations on the lattice, given by:
(7.8) 
where denotes the integration over the direction of the wave vector .
In our implementation of the algorithm we calculate six independent arrays of and use the Fourier transform to obtain the energymomentum tensor . We calculated the projection operators in advance for three families of directions: along edges (), along diagonals of walls () and along the diagonal () of the cubic array of momenta. The integration over directions of the wave vector is substituted by the average of values in distinguished families.
The spectra of GWs obtained from our numerical simulations correspond to a very early time in the evolution of the Universe when the decay of domains walls ended. In order to compare this data with sensitivities of detectors of GWs we must also calculate how GWs emitted from domain walls were stretched by the expansion of the Universe. Solution (7.7) predicts that the energy density of GWs during the radiation domination epoch scales as . Basing on equation (7.2) the ratio of values of the scale factor at the end of the decay of domain walls to the value at the time of equality of matter and radiation energy densities can be expressed by the ratio of proper values of the Hubble constant:
(7.9) 
Values of the Hubble constant at both conformal times can be easily computed. We estimate the value of the Hubble constant at the time of equality of matter and radiation energy densities assuming simple scaling of these densities:
(7.10) 
from present day values of the Hubble constant and fractions of the critical density of matter and radiation . We compute the value of the Hubble constant at the end of the decay of domain walls from parameters of the simulation
(7.11) 
The ratio is relevant for theories in which the radiation domination epoch begins at very low energy scales. In most inflationary scenarios this value is negligible. Finally, substituting equations (7.10) and (7.11) we obtain:
(7.12) 
Furthermore, assuming that the energy density of GWs scales as from the epoch of equality to the present day, we can write:
(7.13) 
where is the present time and is the redshift to the epoch of matterradiation equality. The energy density of gravitational waves is usually presented as a fraction of the critical density .
In order to connect values of the wave vectors on the lattice with the present frequency of GWs we take into account that the wavelength of the GW at any time with the comoving wave vector at the time of the decay of domain walls satisfy
(7.14) 
thus
(7.15) 
where is the present conformal time. Using the previously estimated value of the redshift we get:
(7.16) 
Using this algorithm in lattice simulations raises some complications. The modified eom (A.1) with and cannot be used, because the modification disturbs dynamics of short wavelength fluctuations. For the unmodified eom the width of the domain walls decreases as . The width of domain walls at the end of the simulation expressed in units of the lattice spacing cannot be too small in order to correctly model the profile of walls. On the other hand, the width of domain walls at the initialisation cannot be too large, because many domain walls need to be present on the lattice to properly reproduce the statistical properties of the network. Those requirements significantly restrict the dynamical range of the simulation and consequently only late domain walls () can be reliably investigated in our simulation.
Moreover, as noted in [25], this algorithm produces a spectrum that diverges as (due to the factor in (7.8)) for random initialisation of the field strength. Following [25] we have introduced a cutoff scale of the order of the width of domain walls in the Fourier transform of the initialisation configuration. However, imposing the cutoff decreases the number of degrees of freedom on the lattice (the long wavelength correlations are induced) and the accuracy of the reproduction of the initialisation probability distribution. Thus, the computation of spectra of GWs produced by networks of domain walls from finetuned initial configurations requires accordingly large lattice sizes.
We calculated spectra of GWs produced by decaying networks of Higgs domain walls formed around the conformal time of the order of . We used the value of the suppression scale which leads to the effective potential with nearly degenerate minima. Our simulations were initialised with configurations of the field generated randomly from probability distributions satisfying . We used the value in the initialisation probability distribution. Used values of and correspond to the most stable late networks observed in our scans over generic initialisation conditions.
The figure 11 shows the present spectra of GWs emitted from Higgs domain walls for two values of leading to (red) and (green). Each of these spectra is the average over five simulations performed on the lattice of the size and initialised at with . Simulations ended at the conformal time equal to with . We also assumed that and the value of the Hubble constant at the end of simulations is equal to . The width of domain walls at this time was of the order of . The resulting spectra are peaked at the frequency of the order of and their amplitude is very small.
The figure 12 shows the resulting GW spectra (solid) and sensitivity of present and future detectors LIGO [42], LISA [43] and BBO [44] and the bound on additional energy from CMB/BBN [45, 46] (shaded regions). Spectrum obtained in our previous work [1] for Higgs domain walls in the SM without nonrenormalizable operators is also shown in the figure 12 for comparison. The spectrum of GWs produced by domain walls in the case of nearly degenerate minima has higher amplitude than the spectrum obtained in the case of the pure SM. The difference is a result of longer lifetime of domain walls when nonrenormalizable operators are present. Moreover, the delayed decay of the network effects in the shift of the position of the peak towards lower frequencies. The position of the peak corresponds to the energy scale given by the Hubble radius, thus the later the decay of network of domain walls ends the lower is the frequency of the maximum of the spectrum.
The figure 12 shows that even though the energy density emitted in the scenario with nearly degenerate minima of the effective potential is higher than the previously computed value in the case of pure SM, it is still orders of magnitude lower than predicted sensitivity of planned detectors of GWs. This result excludes the possibility of the detection of GWs emitted by SM domain walls in upcoming years for the generic initial configuration of the Higgs field.
8 Summary
We investigated the evolution of networks of cosmological domain walls of the Higgs field in the presence of corrections from physics beyond the Standard Model. We modelled the influence of heavy states not described by the SM, in the framework of the effective field theory by inclusion of the nonrenormalizable operator suppressed by a scale . We have considered a large range of values of starting from the Planck scale and ending at at which the high field strength minimum turns into a saddle point. We also determined that the electroweak minimum and the high field strength minimum of the RG improved effective potential are degenerate for the value of around .
In section 3 we described properties of domain walls interpolating between the electroweak minimum of the RG improved effective potential and the second minimum located at the high field strengths. We determined the width of domain walls to range from to for the considered range of values of the suppression scale . In addition we calculated the shift of the position of the local maximum separating two minima induced by the inclusion of the operator .
We used numerical lattice simulations based on the PRS algorithm. The same algorithm was used in our previous work [1] and based on that experience we have run simulations for two different values of the initialisation conformal time: and .
Networks of cosmological domain walls could be produced in the early Universe by a variety of processes. In this paper we concentrated on the possibility that Higgs field fluctuations interpolating through the potential barrier separating two minima of the potential were generated during cosmological inflation. Inflationary models predict configurations of the Higgs field well approximated by the Gauss distribution with the standard deviation proportional to the Hubble constant during inflation . However, the mean value of the Higgs field does not change significantly during inflation so it can not be determined just by an inflationary model. In this study we simply assumed the most promising, from the point of view of phenomenology, scenario in which the mean value of the probability distribution vanishes. Our numerical simulations prove that configurations with the initial standard deviation lower than decay to the electroweak vacuum observed today.
Our recent results show that the influence of the nonrenormalizable operator on the evolution of Higgs domain walls was negligible as long as the suppression scale was greater than . For higher our previous results obtained in the pure SM [1] hold. For lower values of dynamics of Higgs domain walls change. With the scale decreasing down to the value at which two minima are degenerate the maximal allowed (by the requirement of the proper final state) standard deviation increases due to the shift of the position of the local maximum to higher field strengths. However, the maximal (allowed by the experimental data) ratio stays constant and approximately equal to . Furthermore, late domain walls i.e. formed in simulations initialised at have longer decay times than early ones (from simulations initialised at earlier conformal times). For values of for which the electroweak minimum is the global minimum of the RG improved effective potential, all simulations with at the initialisation end with the electroweak vacuum as the final state.
We find that networks of Higgs domain walls produced by generic initialisation conditions are highly unstable. The lifetimes of such generic networks decaying into the electroweak vacuum predicted by lattice simulations are shorter than . Moreover, we find no significant increase in the decay time of networks of Higgs domain walls for the effective potential with nearly degenerate minima.
Such short lifetimes call into question the results of previous studies of networks of domain walls with potentials with nearly degenerate minima in which these networks were found to be metastable. However, majority of previous studies considered approximate global symmetries, leading to potentials symmetric up to small corrections. This is not the case in the SM where the scalar potential is highly asymmetric. Thus, our numerical simulations show that previous results can not be straightforwardly extrapolated to the case of the SM with nonrenormalizable operators. The example of the SM proves that degeneracy of minima of the potential does not necessarily implicate metastability of networks of domain walls. Dynamics of domain walls driven by asymmetric potentials may differ significantly from the one observed with nearly symmetric ones.
Our numerical simulations show that lifetimes of networks could be longer if their initial properties were more specific. We considered networks produced by probability distributions satisfying for and . In these finetuned scenarios we observed lifetimes as long as . Initial conditions need to be finetuned in order to compensate for asymmetry of the Higgs effective potential around the local maximum.
Our simulations with finetuned conditions prove the possibility of formation of metastable networks of cosmological domain walls when the minima of the potential are nearly degenerate. Moreover, we showed that these longlived networks of domain walls evolved in the scaling regime satisfying main assumption of semianalytical methods used in [24]. However, semianalytical methods do not accommodate for effects of asymmetry of the potential, other than the difference of values between minima. On the other hand, requirement of removing unphysical, short wavelengths modes from initial configurations limits the usefulness of lattice simulations for computing the energy spectrum of GWs. The calculation for highly finetuned initial conditions requires very large lattice size and we were not able to calculate explicitly the spectrum in our lattice simulations in this case
We obtained longlived networks only for very specific values of the suppression scale , thus both the coupling constant and the initial configuration of the field need to be finetuned. Requirement of the simultaneous finetuning of two, a priori not related, values makes the realisation of this scenario in the nature highly unlikely.
Finally, we computed the energy spectrum of gravitational waves produced during the decay of networks of Higgs domain walls. The longlived networks were of main interest because the later the process of the decay of domain walls ended the larger the present energy density of produced gravitational waves is predicted to be. We used the value of the suppression scale which give the RG improved Higgs effective potential with nearly degenerate minima.
The spectrum of GWs produced by domain walls in case of nearly degenerate minima has higher amplitude than the one previously obtained in the case of the SM without nonrenromalizable operators [1]. The difference is the result of longer lifetime of domain walls when nonrenormalizable operators are present. Moreover, the delayed decay of the network effects in the shift of the peak towards lower frequencies. Its position corresponds to the Hubble radius, thus the later decay of the network ended, the lower is the frequency of the maximum of the spectrum. Even though, the inclusion of nonrenormalizable operators may lead to enhanced energy density of GWs, the one computed for generic initial configurations is still orders of magnitude to low to be detected in the planned detectors.
However, the amplitude of the spectrum of GWs can be larger if the evolution of the network of domain walls took place in the specific conditions. As we pointed out previously in [1], the present energy density of GWs produced by domain walls could be greater if the evolution of the network of domain walls did not take place deep in the radiation era. Models predicting very low scale of the inflation [47] or including new components of energy density which shorten the radiation domination period [48, 49] would result in larger (and lower ). However, this possibility is hard to investigate using lattice simulations due to their small dynamical range and requires semianalytical extrapolations.
Acknowledgments
This work has been supported by the Polish NCN grants DEC2012/04/A/ST2/00099, 2014/13/N/ST2/02712 and 2016/23/N/ST2/03111 and MNiSW grant IP2015 043174 and by the ARC Centre of Excellence for Particle Physics at the Terascale (CoEPP) (CE110001104) and the Centre for the Subatomic Structure of Matter (CSSM).
Appendix A Discretisation of equations of motion
In recent paper we used a discretisation scheme at first proposed in [22]. It was widely used in the past for numerical simulations of dynamics of domain walls for example in [23, 50]. Let us consider the following generalisation of eq. (3.3)
(A.1) 
proposed in [22]. The equation (A.1) with values and combined with used discretisation scheme, proposed in [22], are known as PRS algorithm.
Moreover we have used the adaptive time step method which has been proven in the past to be very useful in cosmological lattice simulations [23, 1]. More specifically, in our simulations the time step is calculated from condition that the maximal value (over the lattice) of the relative error generated in the integration step which can be estimated as:
(A.2) 
must be smaller than some small constant number . Using this method decreases the number of needed time steps leading to reduction of both: a computational time and errors coming from finite accuracy of the representation of floating point numbers.
References
 [1] T. Krajewski, Z. Lalak, M. Lewicki and P. Olszewski, Domain walls and gravitational waves in the Standard Model, JCAP 1612 (2016) 036, [1608.05719].
 [2] ATLAS Collaboration collaboration, G. Aad et al., Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys.Lett. B716 (2012) 1–29, [1207.7214].
 [3] CMS Collaboration collaboration, S. Chatrchyan et al., Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys.Lett. B716 (2012) 30–61, [1207.7235].
 [4] A. Andreassen, W. Frost and M. D. Schwartz, Scale Invariant Instantons and the Complete Lifetime of the Standard Model, 1707.08124.
 [5] A. Salvio, A. Strumia, N. Tetradis and A. Urbano, On gravitational and thermal corrections to vacuum decay, JHEP 09 (2016) 054, [1608.02555].
 [6] Z. Lalak, M. Lewicki and P. Olszewski, Higherorder scalar interactions and SM vacuum stability, JHEP 05 (2014) 119, [1402.3826].
 [7] A. Andreassen, W. Frost and M. D. Schwartz, Consistent Use of the Standard Model Effective Potential, Phys. Rev. Lett. 113 (2014) 241801, [1408.0292].
 [8] D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio et al., Investigating the nearcriticality of the Higgs boson, JHEP 12 (2013) 089, [1307.3536].
 [9] G. Degrassi, S. Di Vita, J. EliasMiro, J. R. Espinosa, G. F. Giudice, G. Isidori et al., Higgs mass and vacuum stability in the Standard Model at NNLO, JHEP 08 (2012) 098, [1205.6497].
 [10] J. Ellis, J. R. Espinosa, G. F. Giudice, A. Hoecker and A. Riotto, The Probable Fate of the Standard Model, Phys. Lett. B679 (2009) 369–375, [0906.0954].
 [11] J. R. Espinosa, G. F. Giudice and A. Riotto, Cosmological implications of the Higgs mass measurement, JCAP 0805 (2008) 002, [0710.2484].
 [12] J. A. Casas, V. Di Clemente and M. Quiros, The Standard model instability and the scale of new physics, Nucl. Phys. B581 (2000) 61–72, [hepph/0002205].
 [13] J. A. Casas, J. R. Espinosa and M. Quiros, Standard model stability bounds for new physics within LHC reach, Phys. Lett. B382 (1996) 374–382, [hepph/9603227].
 [14] J. A. Casas, J. R. Espinosa and M. Quiros, Improved Higgs mass stability bound in the standard model and implications for supersymmetry, Phys. Lett. B342 (1995) 171–179, [hepph/9409458].
 [15] M. Sher, Electroweak Higgs Potentials and Vacuum Stability, Phys. Rept. 179 (1989) 273–418.
 [16] A. Kobakhidze and A. SpencerSmith, The Higgs vacuum is unstable, 1404.4709.
 [17] W. E. East, J. Kearney, B. Shakya, H. Yoo and K. M. Zurek, Spacetime Dynamics of a Higgs Vacuum Instability During Inflation, Phys. Rev. D95 (2017) 023526, [1607.00381].
 [18] A. Joti, A. Katsis, D. Loupas, A. Salvio, A. Strumia, N. Tetradis et al., (Higgs) vacuum decay during inflation, JHEP 07 (2017) 058, [1706.00792].
 [19] A. Hook, J. Kearney, B. Shakya and K. M. Zurek, Probable or Improbable Universe? Correlating Electroweak Vacuum Instability with the Scale of Inflation, JHEP 01 (2015) 061, [1404.5953].
 [20] J. R. Espinosa, G. F. Giudice, E. Morgante, A. Riotto, L. Senatore, A. Strumia et al., The cosmological Higgstory of the vacuum instability, JHEP 09 (2015) 174, [1505.04825].
 [21] L. DarmÃ©, J. Jaeckel and M. Lewicki, Towards the fate of the oscillating false vacuum, 1704.06445.
 [22] W. H. Press, B. S. Ryden and D. N. Spergel, Dynamical Evolution of Domain Walls in an Expanding Universe, Astrophys. J. 347 (1989) 590–604.
 [23] Z. Lalak, S. Lola and P. Magnowski, Dynamics of domain walls for split and runaway potentials, Phys. Rev. D78 (2008) 085020, [0710.1233].
 [24] N. Kitajima and F. Takahashi, Gravitational waves from Higgs domain walls, Phys. Lett. B745 (2015) 112–117, [1502.03725].
 [25] T. Hiramatsu, M. Kawasaki and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, JCAP 1402 (2014) 031, [1309.5001].
 [26] W. Buchmuller and D. Wyler, Effective Lagrangian Analysis of New Interactions and Flavor Conservation, Nucl. Phys. B268 (1986) 621–653.
 [27] B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, DimensionSix Terms in the Standard Model Lagrangian, JHEP 10 (2010) 085, [1008.4884].
 [28] B. Grzadkowski and J. Wudka, Bounds on the HiggsBoson mass in the presence of nonstandard interactions, Phys. Rev. Lett. 88 (2002) 041802, [hepph/0106233].
 [29] C. P. Burgess, V. Di Clemente and J. R. Espinosa, Effective operators and vacuum instability as heralds of new physics, JHEP 01 (2002) 041, [hepph/0201160].
 [30] V. Branchina, E. Messina and M. Sher, Lifetime of the electroweak vacuum and sensitivity to Planck scale physics, Phys. Rev. D91 (2015) 013003, [1408.5302].
 [31] D. Coulson, Z. Lalak and B. A. Ovrut, Nonequilibrium phase transitions and domain walls, in Particles, strings and cosmology. Proceedings, 19th Johns Hopkins Workshop and 5th PASCOS Interdisciplinary Symposium, Baltimore, USA, March 2225, 1995, pp. 429–444, 1995. hepph/9508226.
 [32] Z. Lalak, Towards a solution of the cosmological domain walls problem, in High energy physics, 1996. hepph/9702405.
 [33] J. C. R. E. Oliveira, C. J. A. P. Martins and P. P. Avelino, The Cosmological evolution of domain wall networks, Phys. Rev. D71 (2005) 083509, [hepph/0410356].
 [34] M. Kawasaki and K. Saikawa, Study of gravitational radiation from cosmic domain walls, JCAP 1109 (2011) 008, [1102.5628].
 [35] A. M. M. Leite and C. J. A. P. Martins, Scaling Properties of Domain Wall Networks, Phys. Rev. D84 (2011) 103523, [1110.3486].
 [36] T. Garagounis and M. Hindmarsh, Scaling in numerical simulations of domain walls, Phys. Rev. D68 (2003) 103506, [hepph/0212359].
 [37] P. P. Avelino, J. C. R. E. Oliveira and C. J. A. P. Martins, Understanding domain wall network evolution, Phys. Lett. B610 (2005) 1–8, [hepth/0503226].
 [38] P. P. Avelino, C. J. A. P. Martins and J. C. R. E. Oliveira, Onescale model for domain wall network evolution, Phys. Rev. D72 (2005) 083506, [hepph/0507272].
 [39] A. M. M. Leite, C. J. A. P. Martins and E. P. S. Shellard, Accurate Calibration of the Velocitydependent Onescale Model for Domain Walls, Phys. Lett. B718 (2013) 740–744, [1206.6043].
 [40] C. J. A. P. Martins, I. Yu. Rybak, A. Avgoustidis and E. P. S. Shellard, Extending the velocitydependent onescale model for domain walls, Phys. Rev. D93 (2016) 043534, [1602.01322].
 [41] J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman and J.P. Uzan, Theory and Numerics of Gravitational Waves from Preheating after Inflation, Phys. Rev. D76 (2007) 123517, [0707.0875].
 [42] LIGO Scientific collaboration, J. Aasi et al., Advanced LIGO, Class. Quant. Grav. 32 (2015) 074001, [1411.4547].
 [43] N. Bartolo et al., Science with the spacebased interferometer LISA. IV: Probing inflation with gravitational waves, JCAP 1612 (2016) 026, [1610.06481].
 [44] K. Yagi and N. Seto, Detector configuration of DECIGO/BBO and identification of cosmological neutronstar binaries, Phys. Rev. D83 (2011) 044011, [1101.3940].
 [45] S. HenrotVersille et al., Improved constraint on the primordial gravitationalwave density using recent cosmological data and its impact on cosmic string models, Class. Quant. Grav. 32 (2015) 045003, [1408.5299].
 [46] T. L. Smith, E. Pierpaoli and M. Kamionkowski, A new cosmic microwave background constraint to primordial gravitational waves, Phys. Rev. Lett. 97 (2006) 021301, [astroph/0603144].
 [47] M. Artymowski, Z. Lalak and M. Lewicki, Multiphase induced inflation in theories with nonminimal coupling to gravity, 1607.01803.
 [48] M. Lewicki, T. RindlerDaller and J. D. Wells, Enabling Electroweak Baryogenesis through Dark Matter, JHEP 06 (2016) 055, [1601.01681].
 [49] F. P. Huang, Y. Wan, D.G. Wang, Y.F. Cai and X. Zhang, Hearing the echoes of electroweak baryogenesis with gravitational wave detectors, Phys. Rev. D94 (2016) 041702, [1601.01640].
 [50] A. Lazanu, C. J. A. P. Martins and E. P. S. Shellard, Contribution of domain wall networks to the CMB power spectrum, Phys. Lett. B747 (2015) 426–432, [1505.03673].