Gravitational waves from domain walls in the next-to-minimal supersymmetric standard model

The next-to-minimal supersymmetric standard model predicts the formation of domain walls due to the spontaneous breaking of the discrete $Z_3$-symmetry at the electroweak phase transition, and they collapse before the epoch of big bang nucleosynthesis if there exists a small bias term in the potential which explicitly breaks the discrete symmetry. Signatures of gravitational waves produced from these unstable domain walls are estimated and their parameter dependence is investigated. It is shown that the amplitude of gravitational waves becomes generically large in the decoupling limit, and that their frequency is low enough to be probed in future pulsar timing observations.


Introduction
Supersymmetry (SUSY) is a well-motivated extension of the standard model (SM) of particle physics involving a solution to the hierarchy problem [see, e.g., [1,2] for reviews]. The simplest supersymmetric extension of the SM is called the minimal supersymmetric standard model (MSSM), which consists of only SM particles along with a pair of Higgs doublets and their superpartners. Such a model is also motivated by the unification of the gauge forces and the existence of dark matter candidates. These interesting features will be investigated in various upcoming observations including high-energy physics, astrophysics, and cosmology.
Aside from the rich phenomenology encountered in the MSSM, it suffers from the socalled µ-problem [3], which is originated from the fact that the superpotential of the MSSM contains a term of the form µH u H d (called the "µ-term") at the tree level. Here H u and H d are two Higgs doublet superfields, and µ is some dimensionful parameter. Since the MSSM might be a low energy effective theory derived from a more fundamental theory at higher energies, we naively expect that the magnitude of the µ parameter is of the order of some cutoff scale such as the Planck scale. On the other hand, µ should be of the order of the soft SUSY breaking scale in order to induce the electroweak symmetry breaking appropriately. These considerations seem to revive the hierarchy problem in the MSSM.
The next-to-minimal supersymmetric standard model (NMSSM) [for reviews, see [4,5]] is one of possible solutions to the µ-problem, which introduces an additional singlet superfield S and replaces the µ-term with the coupling of the form λSH u H d . Owing to this coupling, an appropriate magnitude of the µ parameter is dynamically generated when the scalar component of the singlet superfield acquires a vacuum expectation value (VEV) at the electroweak phase transition. However, the introduction of the singlet superfield allows the existence of other dimensionful parameters at the tree level. The simplest way to prohibit such additional dimensionful parameters is to impose a discrete Z 3 -symmetry, and such a model is called the Z 3 -invariant NMSSM, or often abridged as the NMSSM.
The discrete Z 3 symmetry of the NMSSM is spontaneously broken when the scalar components of the Higgs fields as well as the singlet scalar acquire VEVs at the electroweak phase transition, which leads to the formation of domain walls in the universe. The existence of such domain walls turns out to be problematic in the standard cosmology [6], and hence it was argued that the NMSSM is unfavorable in the cosmological point of view [7]. Later, it was pointed out that there are some possibilities to arrange the model to have a small explicit symmetry breaking term, which causes a late time annihilation of the walls [8,9]. Taking account of these possibilities, we deduce that unstable long-lived domain walls can exist in the early epoch of the universe.
In this paper, we point out that a stochastic background of gravitational waves with observably large amplitude can be produced from domain walls in the NMSSM. The amplitude of gravitational waves is determined once we specify the surface mass density of domain walls and their decay time [10,11]. Here we estimate the surface mass density by solving nonlinear field equations for the Higgs sector in the NMSSM to find a planar domain wall solution. On the other hand, the peak frequency of gravitational waves solely depends on the decay time scale. We see that the typical frequency of the gravitational waves corresponds to the band which can be probed via pulsar timing observations. Another interesting finding is that the amplitude of gravitational waves becomes generically large in the so-called decoupling limit, where the couplings between the singlet field and Higgs fields become extremely small. Note that, due to the very weak coupling between the singlet-like states and the MSSM sector in the decoupling regime, the signatures such as those from the collider tend to reduce to those in the MSSM, and we emphasize a complementary role of the gravitational waves as a possible probe to distinguish the NMSSM from the MSSM.
Signatures of gravitational waves in the NMSSM were also investigated in Ref. [12] in the context of the first order electroweak phase transition. Differently from that approach, in this paper we focus on the production of gravitational waves from domain walls, which does not depend on the order of the phase transition. The crucial assumption is the occurrence of the spontaneous breaking of a discrete symmetry followed by the late time collapse of domain walls, which is a definite consequence of the NMSSM with the (approximate) Z 3 -symmetry.
The organization of this paper is as follows. In Sec. 2, we describe the scalar potential in the Z 3 -invariant NMSSM, especially focusing on its Higgs sector. In Sec. 3, we introduce a numerical method to analyze the structure of domain walls, and discuss their cosmological evolution. Then, in Sec. 4, we estimate the amplitude and frequency of gravitational waves produced from domain walls and discuss prospects for future observations. Finally, Sec. 5 is devoted to the conclusions.

Higgs potential in the NMSSM
Let us consider the supersymmetric standard model with an additional gauge singlet superfield S. As was described in Sec. 1, we impose a discrete Z 3 symmetry, under which every chiral supermultiplet Φ transform as Φ → e 2πi/3 Φ, (2.1) in order to forbid the presence of dimensionful parameters. The renormalizable superpotential of the NMSSM allowed by the Z 3 symmetry and the gauge invariance can be written as where λ and κ are dimensionless couplings, and W Yukawa corresponds to the usual Yukawa couplings of the quark and lepton superfields of the MSSM. In the above equations, contraction over gauge indices is understood. For instance, we have for SU (2) L doublet fields, which leads Here and hereafter we use the same symbol for the superfield and its scalar component. From F-term contributions of Eq. (2.2) together with D-term contributions and soft SUSY breaking contributions, we obtain Higgs potential at the tree level: where g 1 and g 2 are the U (1) Y and SU (2) L gauge couplings, respectively, and m 2 Hu , m 2 H d , m 2 S , A λ , and A κ are dimensionful soft parameters. We can rotate away the component H + u by using SU (2) L transformation, and hereafter we take H + u = 0. Furthermore, it is necessary to take H − d = 0 in order to guarantee that the squared masses of charged components become positive. Therefore, in the following we only consider the potential for the neutral components, and simply denote them as H u = H 0 u and H d = H 0 d . The scalar potential becomes where for the Z boson mass m Z ≃ 91.2GeV and the Higgs VEV v ≃ 246GeV. Here we assume that there is no explicit CP violation, and take λ, κ, A λ , and A κ to be real. It is possible to take the VEV of H u to be real and positive by the use of U (1) Y transformation. In general, H d and S are complex, but such an extremum becomes always a local maximum if there is no explicit CP violation [13]. Therefore, we can also take H d and S to be real up to the Z 3 transformations. 1 Then, we parameterize three VEVs by using three real quantities: 1 If we parameterize general three complex VEVs as The no-go theorem of Ref. [13] ensures that these two phases take trivial values at the minimum of the potential. If we fix the U (1)Y gauge such that vu > 0 and φu = 0, this fact implies that there are two more vacua with (φ d , φs) = (2π/3, 4π/3) and (4π/3, 2π/3) in addition to that given by Eq. (2.10). Note that it is also possible to choose the U (1)Y gauge such that these three degenerate vacua are represented by φi = 0, 2π/3, and 4π/3 (i = u, d, s), respectively.
with v u > 0. Furthermore, we see that the potential possesses the following symmetries: λ, [14]. Using these symmetries, we fix the signs of λ and v d to be positive. Other parameters κ, v s , A λ , and A κ can take both signs. The stationary conditions of the potential ∂V /∂H u = ∂V /∂H d = ∂V /∂S = 0 lead to Once we specify the values of v u , v d , and v s , the soft masses m 2 Hu , m 2 H d , and m 2 S are determined from Eqs. (2.11)-(2.13): Hereafter, we use the following six quantities as free parameters of the model: The definition of the µ-parameter in Eq. (2.17) implies that the value of v s becomes much larger than the electroweak scale for λ ≪ 1 even though µ takes a value of O(100)GeV. Such a case is called the decoupling limit, since the coupling between particles associated with the singlet superfield and the MSSM sector becomes negligible for λ → 0. In this limit, dominant terms in the potential (2.6) are given by from which we estimate the value of v s at the minimum as Note that the condition A 2 κ 9m 2 S should be satisfied in order that the extremum given by Eq. (2.20) becomes deeper than the other possible minimum with v s = 0 [15]. Assuming that the magnitudes of |A κ | and |m S | are comparable to the soft SUSY breaking mass scale, from Eq. (2.20) we see that a large value of v s is compatible with κ ≪ 1. Furthermore, it can be shown that κ and λ should be comparable such that the potential does not have unrealistic minima [15,16]. Therefore, the decoupling limit is given by with other dimensionful parameters being fixed.

Domain walls in the NMSSM
. Domain walls are located around boundaries of these three vacua.
In order to find the domain wall solution [i.e., find the spatial configuration of the scalar fields (S(x), H u (x), H d (x)) around the core of the domain wall], we must solve nonlinear field equations for three complex scalar fields S, H u and H d with the potential given by Eq. (2.6). It is straightforward to solve such a system numerically [7,17]: A planar domain wall solution (S(z), H u (z), H d (z)) perpendicular to the z-axis can be found from the following equations Here we use the globally convergent Newton method [18] to solve this boundary value problem. Namely, we first specify the "initial guess" for (S(z), H u (z), H d (z)) within some finite interval −L ≤ z ≤ L, and deform it such that Eq. (3.1) is satisfied with an appropriate boundary conditions specified at z = ±L. Then, we iterate this procedure until Eq. (3.1) is satisfied with an accuracy of O(10 −5 ). For an appropriate choice of the initial guess of the functions (S(z), H u (z), H d (z)), this algorithm converges after a few iteration steps. Figure 1 shows a typical configuration of the domain wall solution. From Fig. 1 (a) we see that the phases of scalar fields continuously change around the core of the domain walls centered at z = 0. We also compute the spatial distribution of the energy density for this solution where we subtracted a constant term such that ρ wall → 0 is satisfied for z → ±∞. As shown in Fig. 1 (b), the energy of the scalar fields is concentrated around the core of the domain wall. Finally, we estimate the surface mass density σ wall of the domain wall by integrating this energy density along the z-axis: For a choice of parameters used in Fig. 1, we obtain σ wall ≃ 1.9 × 10 9 GeV 3 . Although the precise value of σ wall can only be calculated by using the numerical procedure, we can roughly estimate the parameter dependence of σ wall in the simple decoupling limit given by Eq. (2.21). Regarding the fact that the phase of the S field continuously changes across the core of the wall and that the potential can be approximated as Eq. (2.19) in the decoupling limit, we expect that the typical length scale of the spatial variation of the field S around the core of the wall is given by ∆z ∼ |κA κ v s | −1/2 . Furthermore, the estimation shown in Eq. (2.20) implies that all three terms in the right-hand side of Eq. (2.19) are comparable, and hence we roughly estimate the height of the potential around the core of the wall as V ∼ |κA κ v 3 s |. These observations lead to the following scaling behavior of the surface mass density of domain walls: In fact, we confirmed that this behavior holds very accurately for λ, κ < O(10 −4 ) in our numerical results. Eq. (3.5) indicates that the surface mass density of domain walls becomes enormously large in the decoupling limit. We will discuss the observational consequence of this property in Sec. 4.

Domain wall problem and its solution
Let us turn our attention to the cosmological evolution of domain walls. After the formation of domain walls, two kinds of forces act on them. One is the tension force, which smooths out small scale structures to straighten the wall, and the other is the friction force caused by the interactions between particles in the thermal plasma and Higgs fields constituting the wall. The friction force is efficiently induced by particles whose masses are comparable to the electroweak scale, since their interactions with the Higgs fields are relatively strong. However, the number density of these particles drops exponentially as the temperature of the plasma T decreases, and eventually the friction force becomes negligible. The friction effect for the NMSSM domain walls was quantitatively estimated in Ref. [7], and it was shown that this effect becomes irrelevant for T < O(0.1-1)GeV.
Once the friction force becomes negligible, the dynamics of domain walls is mostly determined by the tension force, which acts like p T ∼ σ wall /R wall for a patch of the walls whose typical curvature radius is given by R wall . In this case, it is known that the evolution of domain walls is described by the so-called scaling solution, in which the typical length scales of the system are given by the Hubble radius. In particular, the energy density of domain walls in the scaling regime is given by [10] where A is a dimensionless coefficient called the area parameter. By comparing a naive estimation ρ wall ∼ σ wall /R wall with Eq. (3.6), we estimate the typical curvature radius of the walls in the scaling regime as R wall ∼ t/A, and hence the tension force is given by p T ∼ Aσ wall /t. In the literature, the scaling behavior of domain walls [Eq. (3.6)] is confirmed both numerically and analytically [10,19]. The result of the numerical simulation for the model of a real scalar field with Z 2 symmetry indicates A ≃ 0.8 ± 0.1 [10], while other simulations for the model with N degenerate vacua show that the value of A increases proportionally with N [20]. Since we are interested in the model with three degenerate vacua (i.e. Z 3 invariant NMSSM), we simply estimate the value of A as A ≃ 0.8 × (3/2) = 1.2.
If the domain walls evolve according to the scaling solution, the decrement of the energy density ρ wall ∝ 1/t is slower than that of dusts ∝ a(t) −3 and radiations ∝ a(t) −4 , where a(t) is the scale factor of the universe. This fact implies that domain walls are likely to dominate the energy density of the universe at later times. Such an existence of domain walls leads to many problems in the standard cosmological scenario [6].
A simple solution to the domain wall problem is to introduce an additional term ∆V which explicitly breaks the discrete symmetry in the potential [6,21]. Let us call this additional term a bias. We note that the magnitude of the bias term should be much less than that of the potential around the core of domain walls (∆V ≪ V ) such that the discrete Z 3 -symmetry holds approximately. If such a bias term exists in the potential, it breaks the degeneracy between different vacua, and the energy difference between neighboring domains causes the volume pressure p V ∼ ∆V acting on the walls. Domain walls start to collapse when this pressure p V exceeds the tension force p T ∼ Aσ wall /t, and we can estimate their decay time t dec from the condition p V ≃ p T : If the decay of domain walls occurs before they overclose the universe, we have ρ wall (t dec ) < ρ crit (t dec ) = 3H 2 (t dec )/8πG, where ρ crit (t dec ) is the critical density of the universe at the time t dec , H(t dec ) is the Hubble parameter at that time, and G is the Newton's constant. From this requirement, we obtain Even if domain walls decay before they overclose the universe, the consideration of big bang nucleosynthesis (BBN) leads to another constraint if they are sufficiently long-lived. For domain walls existing around the epoch of BBN, the ratio between their energy density and the entropy density s is given by It is expected that particles radiated from the walls are likely to destroy light elements created during the epoch of BBN. For the abundance of radiated particles, which is comparable with Eq. (3.9), we require t dec 0.01sec in order to avoid the constraints on hadronic energy injection in the standard BBN scenario [22]. This condition gives a lower limit on the magnitude of the bias term: Several remarks on the origin of the bias term [Eq. (3.10)] were made in the literature. A naive expectation is that the discrete symmetry is slightly broken due to gravity, and that a Planck-suppressed higher dimensional operator is responsible for the bias term [23]. However, it turned out that this possibility is unfavorable in the NMSSM [7], since such a Planck-suppressed interaction radiatively induces a tadpole operator [24], which revives the hierarchy problem because of the large VEV of the singlet field. A solution to this problem was proposed in Ref. [8]: It was pointed out that a small tadpole term of the form ∆V ∼ ξm 3 soft S + h.c. (3.11) can be generated by imposing various additional symmetries to constrain the form of nonrenormalizable interactions. Here, m soft is the soft SUSY breaking mass, and the dimensionless coefficient ξ corresponds to the loop suppression factor. The term given by Eq. (3.11) is small enough not to destabilize the singlet field, and it acts as a bias term to remove domain walls. 2 In addition to the solution described above, there is another possibility, where a coupling between the singlet field and additional vector-like matters which are charged under the SM SU (3) C gauge group or a hidden QCD gauge group is introduced [9]. In the presence of such a coupling, the Z 3 -symmetry becomes anomalous for the QCD, and the bias term which destabilize the domain walls is generated due to the instanton effect [28]. In this case, the magnitude of the bias term is estimated as where Λ corresponds to the QCD scale Λ QCD ≈ O(100)MeV if vector-like exotics are charged under SU (3) C , or a scale at which the hidden gauge force becomes strong if they are charged under the hidden gauge group. Once the magnitude of the bias term is specified as Eq. (3.11) or Eq. (3.12), we can estimate the decay time of domain walls according to Eq. (3.7). Since the origin of the bias term is highly model dependent, we take t dec as a free parameter in the rest of this paper.

Gravitational waves from domain walls
From the discussion in the previous section, we see that domain walls exist after the electroweak phase transition, and that they collapse before the epoch of BBN t dec 0.01sec due to the effect of the bias term. In this section, we estimate gravitational wave signatures produced from these unstable domain walls [10,11,29].
The cosmological evolution of domain wall networks is accompanied by various violent processes of collisions and annihilations of them. During these processes, a part of the energy stored in them is released as gravitational waves, which leads to a stochastic gravitational wave background observed today. The amplitude of gravitational waves at the cosmic time t for a certain frequency f can be characterized by the following quantity where ρ crit (t) and ρ gw (t) are the critical density of the universe at the time t and the energy density of gravitational waves at that time, respectively. The spectrum of gravitational waves from domain walls is estimated in Refs. [10,11] by the use of 3-dimensional lattice simulation of the scalar field. According to the recent detailed study [10], the spectrum of gravitational waves from domain walls has a peak at the frequency corresponding to the Hubble radius at the time t dec . The peak amplitude of the gravitational waves at that time is given by [10] Ω gw (t dec ) peak = 8πǫ gw G 2 A 2 σ 2 wall 3H 2 (t dec ) , whereǫ gw is an efficiency parameter of the gravitational radiation. The result of the numerical simulation givesǫ gw ≃ 0.7 ± 0.4 [10]. Assuming that the gravitational radiation is terminated in the radiation dominated era, we estimate the amplitude of the gravitational waves at the present time t 0 as where Ω R h 2 = 4.15 × 10 −5 is the density parameter of radiations at the present time, h is the present value of the Hubble parameter in the unit of 100km·sec −1 Mpc −1 , g S0 = 3.91 is the effective degrees of freedom for the entropy density at the present time, and g 0 = 3.36 and g * are the number of radiation degrees of freedom at the present time and t dec , respectively. Substituting Eqs. (4.4) The peak frequency is given by the Hubble parameter at the decay time: We see that Ω gw h 2 (t 0 ) peak and f (t 0 ) peak are determined by two parameters: σ wall and ∆V [or those appearing in the explicit from of the bias term, i.e., Eq. As was noted in Sec. 3.2, here we use t dec as a free parameter instead of ∆V . In order to estimate the peak amplitude for a given set of the NMSSM parameters [Eq. (2.18)], we must evaluate the value of the surface mass density of domain walls σ wall appearing in Eq. (4.6). Here, we use the numerical method introduced in Sec. 3.1 to estimate σ wall . Figure 2 shows the estimated peak amplitude in two-dimensional parameter space of (λ, κ). In this figure, we also plotted some constraints on the tree-level Higgs potential: One is the condition that the mass matrix of seven Higgs scalars (three CP-even, two CP-odd, and a pair of charged states) do not have any negative eigenvalues, and the other is the requirement that the Higgs potential does not lead to unrealistic minima. We used the analytic expression for unrealistic minima obtained in Ref. [16], and checked whether the depth of the desired minimum with ( S , H u , 2) is deeper than that of possible unrealistic minima for given values of the parameters (λ, κ). 3 The increase in the peak amplitude of gravitational waves for λ ≪ 1 is straightforwardly understood from the relation in the decoupling limit [Eq. (3.5)]. This relation implies that Ω gw h 2 (t 0 ) peak ∝ κ 2 v 6 s t 2 dec ∝ κ 2 λ −6 µ 6 t 2 dec in the limit λ → 0. We see that the decoupling limit is relevant to the observations, since it generically predicts a large amplitude of gravitational waves. Regarding this fact, in Table 1 we choose some bench mark points in the parameter space, whose predictions can be tested in future experimental research.
There are various ongoing and planned gravitational wave experiments, which aim to detect stochastic gravitational wave backgrounds as well as astrophysical sources. Among them, ground-based interferometers such as (Advanced) LIGO [31], VIRGO [32], KAGRA [33], and ET [34] will probe frequency range around f ∼ 100Hz. Space-borne interferometers such as eLISA [35] and DECIGO [36] are also planned to probe lower frequency ranges around 0.1mHz to 10Hz. In addition to these experiments, gravitational wave backgrounds can be probed by using the Pulsar Timing Array (PTA) [37], which is sensitive to the frequency range around 10 −9 -10 −8 Hz. Recently, three projects including PPTA, EPTA, and NANOGrav reported upper limit on the stochastic gravitational wave background of Ω gw h 2 < O(1) × 10 −8 [38].  Figure 2. The peak amplitude of gravitational waves estimated from Eq. (4.6) with the assumption of t dec = 0.01sec in the (λ, κ) plane. Blue solid lines correspond to contours for Ω gw h 2 (t 0 ) peak = 10 −9 , 10 −12 , and 10 −15 . Here we use g * = 10.75, A = 1.2 andǫ gw = 0.7 for numerical coefficients appearing in Eq. (4.6). Other theoretical parameters are fixed as A λ = 150GeV, A κ = −150GeV, tan β = 5 and µ = 200GeV. The dark green region corresponds to the parameters at which the mass matrix of Higgs scalars has some negative eigenvalues, and the light green region corresponds to the parameters at which some unrealistic minimum becomes deeper than the desired minimum. The pink region corresponds to the parameter space where the amplitude of gravitational waves exceeds the present limit Ω gw h 2 < 10 −8 from PTA observations. The gray region corresponds to the parameter space where domain walls dominate the energy density of the universe before t dec [i.e. Eq. (3.8)]. Table 1. Three bench mark points and estimated surface mass density of domain walls, peak amplitude of gravitational waves, and its frequency. Here we used g * = 10.75 for t dec = 10 −2 sec and g * = 68.8 for t dec = 10 −6 sec [30] to estimate Ω gw h 2 (t 0 ) peak .
10 −2 sec 10 −2 sec 10 −6 sec σ wall 1.96 × 10 4 TeV 3 1.96 × 10 2 TeV 3 1.96 × 10 8 TeV 3 Ω gw h 2 (t 0 ) peak 4.66 × 10 −9 4.66 × 10 −13 2.51 × 10 −9 f (t 0 ) peak 1.02 × 10 −9 Hz 1.02 × 10 −9 Hz 1.02 × 10 −7 Hz As shown in Fig. 2, this present upper limit already excludes some parameter region for t dec = 0.01sec. The future PTA observations such as FAST [39] and SKA [40] will further improve the sensitivities. Figure 3 shows the schematics of the spectrum of gravitational waves predicted from three points enumerated in Table 1 and sensitivities of planned experiments. The recent result of numerical simulations [10] indicates that the spectrum of gravitational waves from domain walls scales like Ω gw ∝ f 3 for f < f peak and Ω gw ∝ f −1 for f > f peak , and we also plot this expectation in the figure. We see that the point I leads to the peak at the frequency f ∼ 10 −9 Hz, which can be probed in PTA observations such as SKA. On the other hand, the prediction of the point II indicates that the amplitude gets smaller as the values of λ and κ increase. Interestingly, the spectrum of gravitational waves predicted from the point III might be probed by DECIGO as well as SKA. From these illustrations, we expect that the future PTA observations and space-borne interferometers can probe the signature of gravitational waves from NMSSM domain walls.  [42,43] [see also Ref. [10]]. The sensitivity of SKA is taken from Ref. [44].

Conclusions
In this paper, we estimated the gravitational wave signatures from domain walls created at the electroweak phase transition in the NMSSM. The surface mass density of domain walls σ wall is estimated by solving nonlinear equations for Higgs fields numerically. By using the estimated value of σ wall , the amplitude of relic gravitational waves is calculated via Eq. (4.6) and its parameter dependence is explored. We find that a large amplitude of gravitational waves, which scales as Ω gw h 2 peak ∝ κ 2 v 6 s t 2 dec , is predicted in the decoupling limit. Furthermore, the peak frequency is determined by the decay time of domain walls [Eq. (4.7)], which is related to the magnitude of the bias term ∆V via Eq. (3.7). Such gravitational wave signatures can be probed in the future experiments with improved sensitivities.
We note that the decay time of domain walls t dec is restricted to a certain interval, from the epoch of the electroweak phase transition to that of BBN. The typical frequency of gravitational waves corresponding to such time scales is very low, which is suitable for PTA observations. Since the amplitude of gravitational waves becomes large in the decoupling limit, we can use the results of PTA observations as a lower limit on the couplings (λ and κ) appearing in the Higgs sector of the NMSSM. For instance, the present PTA bound Ω gw h 2 O(10 −8 ) already puts a lower limit λ > 10 −4 -10 −3 for t dec = 0.01sec (see Fig. 2), while this constraint becomes irrelevant if the decay time of domain walls t dec is earlier than O(0.01)sec. Future observations will probe a wider range of the parameters in the weak coupling regime, which can be used to complement the results of other experimental studies of the NMSSM.