Interplay between the muon $g-2$ anomaly and the PTA nHZ gravitational waves from domain walls in next-to minimal supersymmetric standard model

With some explicitly $Z_3$ breaking terms in the NMSSM effective superpotential and scalar potential, domain walls (DWs) from spontaneously breaking of the discrete symmetry in approximate $Z_3$-invariant NMSSM can collapse and lead to observable stochastic gravitational wave (GW) background signals. In the presence of a hidden sector, such terms may originate from the geometric superconformal breaking with holomorphic quadratic correction to frame function when the global scale-invariant superpotential is naturally embedded into the canonical superconformal supergravity models. The smallness of such mass parameters in the NMSSM may be traced back to the original superconformal invariance. Naive estimations indicate that a SUSY explanation to muon $g-2$ anomaly can have tension with the constraints on SUSY by PTA data, because large SUSY contributions to $\Delta a_\mu$ in general needs relatively light superpartners while present $\Omega_{gw}^0$ can set the lower bounds for $m_{soft}$. We calculate numerically the signatures of GW produced from the collapse of DWs and find that the observed nHZ stochastic GW background by NANOGrav, etc., can indeed be explained with proper tiny values of $\chi m_{3/2}\sim 10^{-14}{\rm eV}$ for $\chi S^2$ case (and $\chi m_{3/2}\sim 10^{-10}{\rm eV}$ for $\chi H_u H_d$ case), respectively. Besides, there are still some parameter points, whose GW spectra intersect with the NANOGrav signal region, that can explain the muon $g-2$ anomaly to $1\sigma$ range.

Low-energy SUSY is the most promising new physics beyond the SM of particle physics.The SUSY framework can solve almost all the theoretical and aesthetic problems that bother the SM.In particular, the discovered 125 GeV Higgs scalar lies miraculously in the small '115 − 135' GeV window predicted by the low-energy SUSY, which can be seen as a strong hint of the existence of low-energy SUSY.NMSSM can elegantly solve the µ problem of Minimal Supersymmetric Standard Model (MSSM) with an additional singlet sector.In addition, with additional tree-level contributions or through doublet-singlet mixing, NMSSM can accommodate easily the discovered 125 GeV Higgs boson mass with low electroweak fine tuning.The Z 3 -invariant NMSSM is the most predictive version of NMSSM realizations with less new input parameters.However, the VEVs of Higgs fields S, H u , H d at the electroweak phase transition era will also trigger the breaking of the Z 3 discrete symmetry, potentially generated cosmologically problematic DWs because the DWs energy density decreases slower than radiation and matter, and would soon dominate the total energy density of the universe.As long as the discrete symmetry is exact, such topological configurations are stable.However, it was argued that DWs with a tension as large as σ > O(MeV 3 ) must not exist in the present universe, which is known as Zel'dovich-Kobzarev-Okun bound [107].Fortunately, the degeneracy of the minima from the discrete symmetry can be lifted by an energy bias in the potential, causing the collapse of DWs at the early time of the Universe and produce large amounts of GWs.Such GWs can form a stochastic background that persist to the present Universe, potentially being probed by recent PTA experiments.
There are several terms in the general NMSSM superpotential that can explicitly violate the Z 3 discrete symmetry to trigger the decay of DWs [108,109].We argue that the smallness of the quadratic terms may be traced back to the original superconformal invariance for the matter part of the supergravity action when the global scale-free NMSSM is embedded to supergravity.Unlike the pure supergravity part, which breaks superconformal symmetry after the gauge fixing, the matter part alone remains superconformal.Additional real parts of the holomorphic functions in the frame function are adopted to break the superconformal symmetry, which result in the scale-invariance violation quadratic terms (such as H u H d or S 2 terms) in the superpotential in the presence of a hidden sector.
It is well known that the theoretical prediction of the muon anomalous magnetic moment a µ ≡ (g − 2) µ /2 for SM has subtle deviations from the experimental values.In fact, combining the recent reported FNAL muon g − 2 measurement with the previous BNL+FNAL results [110][111][112], the updated world average experimental value of a µ [113] has a 5.1σ deviation from the SM predictions [114] with This paper is organized as follows.In Sec 2, we discuss the formation of DWs from the spontaneous breaking of Z 3 discrete symmetry.In Sec 3, the collapse of DWs from explicitly Z 3 -violation terms are discussed.Such small Z 3 violation terms can be the consequence of superconformal violation in the matter part of supergravity in the presence of a hidden sector.In Sec 4, the GWs signals related to the collapse of DWs are discussed.The interplay between the muon g − 2 anomaly and the PTA nHZ GWs signals from DWs are also discussed.Sec 5 contains our conclusions.

DWs from Z 3 -invariant NMSSM
As noted previously, NMSSM is well motivated theoretically to solve the µ problem.A bare µ term µH u H d in the NMSSM superpotential is forbidden if we impose a discrete Z 3 symmetry, under which the Higgs superfields S, H u , H d transform nontrivially.The Z 3 invariant superpotential couplings are given by [105,106] with 2) The soft SUSY breaking parameters can be given as by proper SUSY breaking mechanisms, such as AMSB [115] and Yukawa mediation mechanisms [116,117].From the superpotential and the soft SUSY breaking terms, one obtains the Higgs potential where g 1 and g 2 denote the U (1) Y and SU (2) L gauge couplings, respectively.Electroweak symmetry breaking (EWSB) are triggered by the VEVs is generated by the VEV s ≡ ⟨S⟩.Such Higgs VEVs also break spontaneously the Z 3 discrete symmetry and lead to the formation of DWs.
As the Z 3 -NMSSM scalar potential is invariant under the Z 3 transformation, the scalar potential has three degenerate minima that can be parametrized by We use the conclusion in [118] that the true minimum of the potential does not spontaneously break CP as the VEVs can always be made real by an appropriate field redefinition, up to the existence of three degenerate vacua related to each other by Z 3 transformations.
DWs are located around boundaries of these three vacua.A planar domain wall solution perpendicular to the z-axis can be found by solving numerically the equation of motions for S, H u , H d fields with the boundary conditions such that the configuration interpolates smoothly between two vacua at z → ±∞.We adopt linearized interpolation curves for the two asymptotic domains as our initial guess of the solution and iterate these procedures until the solution converges.From the configuration of the DWs, the spatial distribution of the energy density can be calculated by in which the boundary energy constant is subtracted such that ρ w → 0 is satisfied for z → ±∞.The tension of the DWs, which measures the energy stored per unit area on the wall, can be calculated from the energy-momentum tensor of the static solution.The calculation can be reduces to the integration of the energy density along the z-axis We show in Table2 the tension and energy density of DWs from Z 3 -NMSSM for some benchmark points.We know that the validity of perturbation theory up to the GUT scale implies λ < 0.7 ∼ 0.8 [105] at the weak scale.Besides, the EWSB condition requires that (2.9) So, the singlet VEV s always takes values of order the soft SUSY breaking scale, which is always much larger than the EW scale VEV v u and v d and act as the dominant spontaneously Z 3 discrete symmetry breaking source.In this case, as noted in [108], the dominant contribution to the energy density and the tension of DWs comes from the singlet relevant terms in the scalar potential.From the equation of motion for the corresponding field, we can estimate the gradient of the field as So, the characteristic length (thickness) for varying ϕ(z) from two adjacent minima −v/2 to v/2 in NMSSM can be estimated to be while the energy density of the DWs is 12) The tension of DW can be estimated in terms of the height of the potential energy V 0 separating the degenerate minima and the thickness δ by the relation σ ∼ δV 0 , giving σ ∼ (κA k s 3 )/A κ ∼ A 3 κ /κ 2 .Such an estimation can be helpful to understand some of the subsequent numerical results.
When the exponentially damped friction force of DWs becomes irrelevant after the temperature of the background radiations becomes less than the mass of particles that interact with DWs, the dynamics of DWs is dominated by the tension force, which stretches them up to the horizon size.The evolution of DWs in this regime can be described by the scaling solution, which states that their energy density evolves according to the simple scaling law and there is about one DW per Hubble horizon.The energy density of DWs in the scaling regime decreases much more slowly than that of cold matters and radiations so that they gradually dominate the energy density of the Universe, drastically altering the subsequent evolution of the Universe.

Collapse of DWs from explicit Z breaking
The formation of DWs from spontaneously discrete Z 3 breaking can be problematic in cosmology if they persist till now.They need to collapse before primordial Big Bang nucleosynthesis (BBN).If the discrete symmetry is broken explicitly by 1/M P suppressed interactions, such operators can lead to a divergent two-or three-loop diagrams that contribute a term linear in S, destabilizing the Planck/weak hierarchy unless the coupling is tiny.However, such tiny couplings conflict with the previous constraints from nucleosynthesis [119].Therefore, explicitly Z 3 breaking model should have nonvanishing µ 0 ̸ = 0 or µ ′ ̸ = 0 with κ ̸ = 0 in the low-energy effective superpotential and, at the same time, the absence of the 1/M P suppressed S 4 /M P operator1 .Although the superpotential will not introduce new terms by the nonrenormalization theorem, such terms and soft SUSY breaking terms always appear in the low-energy effective action for the light superfields [121] after integrating out the messenger/sequestered sector for the SUSY breaking.However, such µ 0 and µ ′ terms can again lead to the existence of tadpole diagrams for S, possibly causing the shift in the potential for S to slide to large values far above the weak scale.In fact, a more general conclusion in [119] states that, to ensure a model with singlets will be natural, typical symmetry (such as the gauged R-symmetry or modular symmetry for target space duality [122] instead of ordinary gauge symmetry) is needed to forbid odd-dimension terms in the Kahler potential K and even-dimension terms in the superpotential Ŵ , while any extra odd-dimension operators in Ŵ or even-dimension operators in K are not harmful to the gauge hierarchy.So, the low-energy effective µ ′ S 2 and µ 0 H u H d operators should be absent in the original superpotential Ŵ for supergravity to avoid problematic singlet tadpole problems and they can only come from the Kahler potential.Such small µ 0 and µ ′ terms can be the natural consequence of small geometrical breaking of larger local superconformal symmetry from the bilinear dimensionless holomorphic terms in the frame function for the canonical superconformal supergravity (CSS) model [123], when the scale-invariant superpotential is embedded into it.

Explicitly Z 3 violation term from geometrical superconformal breaking
The form of Poincare supergravity can be obtained from the underlying superconformal theory, which in addition to local supersymmetry of a Poincare supergravity-has also extra local symmetries: Weyl symmetry, U (1) R symmetry, special conformal symmetry and special supersymmetry.General theory of supergravity in an arbitrary Jordan frame was derived in [124] from the SU (2, 2|1) superconformal theory after gauge fixing.It was noted in [123] that the embedding of the globally superconformal theory into supergravity in the Jordan frame can be realized by simply adding the action of the global SUSY, interacting with gravity with conformal scalar-curvature coupling, to the action of supergravity.The scalar-gravity part of the N = 1, d = 4 supergravity in a generic Jordan frame with frame function Φ(z, z) and superpotential W (z) is given by [123] with and A µ depends on scalar fields as follows: For a particular choice of the frame function, the kinetic scalar terms are canonical when the on-shell auxiliary axial-vector field A µ vanishes.This requires that These J (z) + J (z) terms in the frame function not only break the continuous R symmetry, but also break the discrete Z 3 symmetry.
In addition to the pure supergravity part in the total supergravity action, which breaks superconformal symmetry, the part of the action describing chiral and vector multiplets coupled to supergravity has a much larger local superconformal symmetry.Such a superconformal invariance needs to be broken to make most of the particles massive.However, this symmetry is broken down to the local Poincare supersymmetry only by the part of the action describing the self-interacting supergravity multiplet, sometimes requiring additional superconformal symmetry breaking sources in the matter sector.It was proposed in [125] that the superconformal symmetry of the matter multiplets in the supergravity action can be broken geometrically without introducing dimensional parameters into the underlying superconformal action, in which the moduli space of chiral fields including the compensator field is no longer flat.
We adopt the setting of Jordan frame NMSSM with the following real function that describes the Kahler manifold of the embedding space, which leads to the superpotential in Jordan frame with G α β the matter part of the inverse metric G I J of the enlarged space including the compensator.The general expression of G α β can be calculated to be [125] for the real function So, for our case, we have The SUSY preserving part of the scalar potential receives an additional term (3.12) We can see that the leading Z 3 breaking term in this scalar potential that proportional to (S 2 + S * 2 ) is suppressed by 1/M 4 P .Similarly, for real function N (X, X) of the form in the Jordan frame for NMSSM, we have So, the scalar potential for NMSSM is given by with the correction Therefore, the leading Z 3 breaking term that proportional to (H u H d + Hu Hd ) in the scalar potential can be seen to be suppressed by 1/M 4 P .The frame function Φ(z a , za ) and consequently the Kahler potential can relate to N (X, X) by the gauge fixing with N (X, X) = Φ(z a , za ) = −3e −K(z,z)/3 .In many popular SUSY breaking mechanisms, the total superpotential can be divided into the hidden sector and visible sector with W = W vis + W hid for ⟨W vis ⟩ ≪ ⟨W hid ⟩, as the VEV is obviously dominated by the fields responsible for SUSY breaking even though all fields in the theory are summed over.So, under Kahler transformation to render the new Kahler potential canonical the superpotential will receive a correction of the from Here the gravitino mass m 3/2 is For J(z) = χS2 (or χH u H d ), the superpotential will receive (non-Planck scale suppressed) Z 3 breaking terms χm 3/2 S 2 (or χm 3/2 H u H d ), respectively.Note that unpreferable tadpole term ξ F S will not be present in such Z 3 explicitly breaking superpotential because no dimensional parameters are introduced in N (X, X), which can be seen as an advantage of this CSS approach.
After the Kahler transformation with new frame function Φ ′ (z, z), the action in the Jordan frame is given by with V ′ J the scalar potential for the superpotential e −J W .It is obvious that the matter part of supergravity action contains nonminimal gravitational couplings of scalars to the Ricci scalar curvature in additional to the bosonic part of the supergravity action.To absorb the nonminimal gravitational couplings, we need to rescale the metric.Under the transformation gµν = e 2σ(x) g µν , (3.22) the Ricci scalar changes as So, the metric can be rescaled as followings to change the Jordan frame into the Einstein frame 2 .After metric rescaling, the scalar potential in the Einstein frame is given by with V ′ J (e −J W ) the SUSY scalar potential for the superpotential e −J W .General discussions for nonminimal gravitational couplings involving the complex scalars can be found in our previous work [126], in which the polar variables R, Θ can be used to simplify the transformations.As anticipated, in small field values ϕ ≪ M P , the new variables in the Einstein frame can reduce to the original field variables.The m 3/2 J(z) term in the superpotential can give the dominant (non-Planck scale suppressed) Z 3 breaking contribution in the Einstein frame scalar potential in comparison to the Z 3 breaking terms in (3.12) and (3.16).

Additional explicitly Z 3 breaking terms after SUSY breaking
In the CSS approach with geometrical superconformal breaking in the matter sector, no tadpole terms in the superpotential and scalar potential will be generated after the Kahler transformation with SUSY breaking.The mediation mechanism of SUSY breaking that is responsible for the soft SUSY breaking parameters of NMSSM fields should also not reintroduce the tadpole terms.Unlike some gravity mediation mechanism, which can potentially trigger the generation of tadpole terms again, the tadpole terms in gauge mediated SUSY breaking (GMSB) mechanism can be forbidden by some discrete symmetry for the messenger sector.For example, the economical realization [127,128] (without the introduction of additional gauge fields) can be adopted to avoid the tadpole problem, which can preserve the Z 3 discrete symmetry in the messenger sector by introducing double families of messengers (in 5 ⊕ 5 representation of SU (5)).In this case, the general soft SUSY breaking parameters in the NMSSM, such as B and B ′ , can be estimated to leading order as for the VEV of spurion ⟨X⟩ = M m + θ 2 F X .For a general messenger sector with proper messenger-matter interactions involving the H u or H d fields, the B parameter will receive similarly contribution of order the soft SUSY breaking masses In our subsequent studies, we will not specify concretely which SUSY breaking mechanism is used.The values B ≃ B ′ ∼ m sof t with m sof t = √ m t1 m t2 are adopted in our numerical results for such explicitly Z 3 breaking soft SUSY parameters.

Collapse of DWs and GWs
The χm 3/2 S 2 related terms in the scalar potential break explicitly the Z 3 symmetry.The scalar potential relevant to such a term is given by with µ ′ = 2χm 3/2 .We will neglect the subleading M P suppressed Z 3 breaking terms in our subsequent studies for the collapsing of DWs.The vacuum energy differences between v i and v j are defined as We can calculate the energy difference between adjacent vacuum Requiring the v 1 vacua to be the true minimum leads to the constraints ∆ 1,ω V ̸ Z 3 ;S 2 < 0 and ∆ 1,ω 2 V ̸ Z 3 ;S 2 < 0, which can be satisfied with The degeneracy of the two minimum is the consequence of CP conjugation of the ω and ω 2 vacua.When the conditions in Eq.(4.4) are satisfied, we have the order of the vacuum energies for the biases.With bias terms, the energy difference between the neighboring vacua acts as a volume pressure p V ∼ V bias on the wall, which tends to shrink the higher-energy domains until the whole space is filled by the only true vacuum, eventually collapsing the wall network.The collapse of DWs happens when this pressure force becomes greater than the tension force p T ∼ σ/R wall .In our case with two degenerate false vacua, the volume pressures from (V bias ) 1,ω 2 and (V bias ) 1,ω are of the same strength, and they lead the true v 1 vacuum to expand against the false v ω and v ω 2 vacua regions, reducing at the same rate the occupancy of such false vacuum states and leading to the fragmentation of the wall network, which will collapse.
If the conditions in Eq.(4.4) are not satisfied, the potential has two degenerate true minima v ω and v ω 2 .To spoil the problematic degeneracy between the ω and ω 2 vacua, one can introduce additional messengers ψ ′ , ψ′ , which are charged under some hidden gauge group SU (N ) h and Z 3 .The Z 3 discrete symmetry can be anomalous under SU (N ) h , leading to additional contributions to bias of order ∆ H V ∼ λ 4 H [129] so as that the degeneracy between the ω and ω 2 vacua is spoiled (assuming that V v ω 2 ≳ V vω ).Then, the vacuum energies and the biases satisfy In such a scenario, due to the fact that the true vacuum v ω is nearly degenerate with the false vacuum v ω 2 , the annihilations driven by (V bias ) 1,ω 2 and (V bias ) 1,ω are of similar strengths, increasing almost identically the v ω and v ω 2 regions by volume pressure to diminish the original false vacuum v 1 region.As v ω 2 is still a false vacuum, the tiny (V bias ) ω,ω 2 driven volume pressure on v ω 2 can slowly push the walls to increase the true vacuum v ω region.So, after the period with quickly increment of the v ω 2 region driven by the (V bias ) 1,ω 2 from original v 1 region, the v ω 2 region tends to diminish slowly by the tiny (V bias ) ω,ω 2 driven volume pressure.In this case the decay time for DWs is dominated by the slow (V bias ) ω,ω 2 driven collapse.Collapse of DWs in such a scenario with vacuum structure satisfying (4.5) is sensitive to the small contributions (V bias ) ω,ω 2 that lift the degeneracy between v ω and v ω 2 .Not specifying the details of the degeneracy lifting mechanism, we will not discuss such a scenario in this work and only concentrate on the scenario with v 1 the true vacuum that was discussed in previous paragraphs.The collapse of DWs happens when the volume pressure force becomes greater than the tension force.Assuming that DWs have reached the scaling regime, we can estimate the decay time with A ∼ 1.2 for N = 3 and C ann a coefficient of O(1) that takes value C ann = 5.02 ± 0.44, based on the simulation of axion models with N = 3.We can estimate the decay time as t ∼ C 0 (χm 3/2 ) −1 by the following observations: with C 0 , C 1 some constant coefficients.So, the decay time can be dominantly determined by the value of χm 3/2 .It is well known that DWs cannot be formed from the beginning if the bias terms are sufficiently large.The prediction of percolation theory gives a necessary condition for the formation of large scale DWs with the presence of bias term [132] We check that this constraint is safely satisfied for the NMSSM scalar potential with small bias terms, for example, for µ ′ S 2 /2 case with the choice χm 3/2 ∼ 10 −12 eV adopted in our numerical studies.Similarly, the µ 0 ≡ χm 3/2 H u H d term that explicitly violates the Z 3 discrete symmetry will also cause the collapse of DWs.(We have µ ef f = λs + µ 0 ≈ λs as the term λs gives the dominant contribution to µ ef f .)The scalar potential relevant to such a term is given by We can calculate the energy difference between adjacent vacuum From the electroweak symmetry breaking condition in Eq.(2.9), we have µ ef f ≫ m sof t tan β so that the contribution to ∆ 1,ω;ω 2 V from soft SUSY breaking B term is always negligible.Requiring the v 1 vacua to be the true minimum leads to the constraints µ 0 < 0. Again, in the scaling regime, the decay time of DWs can be estimated as Similarly, the decay time can be dominantly determined by the value of χm 3/2 .The constraint (4.7) is also safely satisfied for the NMSSM scalar potential with tiny explicit Z 3 breaking coefficient χm 3/2 ∼ 10 −12 eV for µ 0 H u H d case.
The collapse of DWs should occur before they overclose the Universe and their decay products should at the same time not destroy light elements created at the epoch of BBN.So, we have the constraints t ≲ min(0.01sec, t dom ) = min(1.4× 10 13 eV −1 , t dom ) , (4.11) with The BBN constraint is always more stringent than that of t dom unless σ ≳ 10 6 TeV 3 .

GWs from the collapse of DWs
The collisions of the DWs can deviate from the spherical symmetry, so the energy stored inside DWs can not only be released into light degrees of freedom but also into GWs when DWs collapse.The intensity of the GWs decreases as the Universe expands since the amplitude is redshifted for each mode.So, one may use a dimensionless quantity Ω gw (f ) to characterize the intensity of GWs.The density parameter Ω gw (f ) of the GWs can be defined as the fraction between the GWs energy density per logarithmic frequency interval and the total energy density of the Universe with ρ gw as the energy density of the GWs, ρ c the critical energy density and f the frequency of GWs.As the GWs are predominantly sourced by horizon-scale structures in the DW network, the GWs emitted are expected to peak at the frequency corresponding to the horizon scale at the decay time.
Due to the redshift by the cosmic expansion, the present density parameter Ω 0 gw and the frequency of the GWs f 0 can be obtained from the quantities at the formation of the GWs (denoted with ′ * ′ ) as The peak frequency today can be estimated by with g * and g * s the effective relativistic degrees of freedom for the energy density and the entropy density, respectively.Here we use the following formula to transform the temperature into the Hubble parameter in the radiation dominant era with M P = 1.22 × 10 19 GeV.
For s, m sof t ∼ O(1) TeV, λ, κ ∼ O(1) and tan β ∼ O (10), the value of χm 3/2 should be of order 10 −12 eV for σ ∼ O(TeV) 3 .Given that the gravitino mass m 3/2 cannot be much lighter than O(1) eV, the χ parameter should be much smaller than unity.Such tiny value can be the consequence of higher-dimensional operators in the Kahler potential with U (1) R symmetry.For example, assuming R(Φ) = 16/3 and R(Φ ′ ) = 4, the Kahler potential adopted in [133]   can naturally lead to χ ∼ 10 −24 for m 3/2 ∼ O(1)GeV and α, α ′ ∼ O(1) without the tadpole problem, when the R-symmetry breaking scale lies of order ⟨Φ⟩ ∼ ⟨Φ ′ ⟩ ∼ 10 6 GeV.Such Kahler potential can originate from the real function N (X, X) after gauge fixing.The power of the GWs radiation can be estimated as ρ gw ∼ P t/t 3 ∼ GAσ 2 by the quadrupole formula, given by P ∼ G ...
wall /t 2 with M wall = σAt 2 the energy of DWs.The quadrupole formula cannot be directly applied to DWs, since it is only valid in the far-field regime [134].To obtain the peak amplitude of GWs produced by long-lived DWs, we use the fact that numerical simulations show that the value almost keeps the constant value εgw ≃ 0.7 ± 0.4 [135] after DWs enter into the scaling regime.So, the peak amplitude at the annihilation time of DWs can be estimated as Here the production of GWs is assumed to suddenly terminate at H = H * and happens during the radiation dominated era.The present density parameter Ω 0 gw can be estimated as [134] Ω gw (t 0 ) peak h 2 = Ω rad h 2 g * (T ann ) g * 0 g * s0 g * s (T ann ) with Ω rad h 2 = 4.15 × 10 −5 the density parameter of radiations at the present time.The constraint t ≲ t dom can set an upper bound for Ω gw (t 0 ) peak h 2 with for the case with σ ≳ 10 6 TeV 3 .The full GW spectrum today can be estimated by In previous discussions, it is assumed that the annihilation of DWs happens during the radiation dominated era.If it happens before reheating, the above estimations can be modified accordingly.Numerical results indicate that such a possibility requires very low reheating temperature of order MeV to explain the recent nHZ GWs background signals.Low reheating temperature of order GeV can be preferred to avoid the overproduction of very light gravitino and the increase of entropy due to the production of moduli fields [136].However, this low temperature seems to be problematic with baryogenesis, according to which it is commonly assumed that the reheating temperature is at least of the order the electroweak scale (10 2 GeV) [137].Although there are some mechanisms that can explain the baryon asymmetry at very low temperatures, reheating temperature of order MeV still seems not favored.

Tension between nHz PTA GW background data and muon g − 2 anomaly
The PTA observations on the frequency of stochastic GW background can set stringent constraints for the collapse time of DWs, which can be mainly determined by the value of χm 3/2 [see the discussions below Eq.(4.6)].With the very narrow range of the DW collapse time derived from PTA data, the upper bound for the present density parameter Ω 0 gw can set an upper bound for the tension σ of DWs by the formula (4.21), which can be further translated into the upper bound for the soft SUSY breaking parameters m sof t .We can estimate the upper bound of m sof t to be m sof t ≲ 200TeV.Such an upper bound on m sof t can shed light on the UV SUSY breaking mechanisms.
Most importantly, given the BBN upper bounds for the collapse time of DWs, the lower bound for the present density parameter Ω 0 gw can set the lower bounds for m sof t , which can be estimated to be of order TeV.Such lower bounds can have tension with the SUSY explanation of the muon g − 2 anomaly because large SUSY contributions to ∆a µ in general need relatively light superpartners.
We know that the SUSY contributions to muon g − 2 are dominated by the charginosneutrino and the neutralino-smuon loops in MSSM.At the leading order of tan β and m W /m SU SY , they are evaluated as [138] masses and µ the Higgsino mass, respectively.The loop functions are defined as which are monochromatically increasing for x > 0, y > 0 with 0 ≤ f C,N (x, y) ≤ 1.They satisfy f C (1, 1) = 1/2 and f N (1, 1) = 1/6 in the limit of degenerate masses.The SUSY contributions to the muon g − 2 will be enhanced for small soft SUSY breaking masses and large value of tan β.It can be seen from the previous formulas that large SUSY contributions to ∆a µ require light sleptons, light electroweakinos and large tan β.However, current LHC experiments have already set stringent constraints on colored sparticles, such as the 2.2 TeV bound for gluino and the 1.4 TeV bound for squarks, making the explanation of the recent muon g − 2 data fairly nontrivial.SUSY explanations of the muon g − 2 anomaly can be seen in the literatures; see [139][140][141][142][143][144][145][146][147][148][149][150][151][152][153][154] and [155][156][157][158][159]. The inclusion of the singlet component in NMSSM will in general give negligible contributions to ∆a µ because of the suppressed coupling of singlino to the MSSM sector.However, the lightest neutral CP-odd Higgs scalar could give non-negligible contributions to a µ if it is quite light [160].
The positive two-loop contribution is numerically more important for a light CP-odd Higgs at approximately 3 GeV and the sum of both one-loop and two-loop contributions is maximal around m a 1 ∼ 6 GeV.However, such a light a 1 is constrained stringently by various lowenergy experiments, such as CLEO and B-physics.
Given the estimated O(TeV) lower bounds for soft SUSY breaking parameters by PTA GW data, it is interesting to survey numerically if the parameter regions allowed by the muon g − 2 explanation in Z 3 -NMSSM can be consistent with the DW explanation of such stochastic GW background signals.

Numerical results
We calculated numerically the tension of DWs, generated from the spontaneously discrete symmetry breaking in the Z 3 -invariant NMSSM, within the parameters regions allowed by low-energy experimental data, including the following constraints other than those already encoded in the NMSSMTools_5.6.2 package (i) The CP-even component S 2 dominated combination of H u and H d doublets corresponds to the SM Higgs, which should lie in the combined mass range for the Higgs boson, 122GeV < M h < 128GeV [161,162].We adopt the uncertainty ±3 GeV instead of default 2 GeV because large λ may induce additional O(1) GeV correction to M h at the two-loop level [163].
(ii) Bounds for low mass and high mass resonances at LEP, Tevatron, and LHC are taken into account by using the package HiggsBounds-5.5.0 [164]and HiggsSignal-2.3.0 [165].
(v) Vacuum stability bounds on the soft SUSY breaking parameters, including the semianalytic bounds for nonexistence of a deeper charge/color breaking minimum [177] and/or a metastable EW vacuum with a tunneling lifetime longer than the age of the Universe [178].
(vi) The relic abundance of the dark matter (DM) should be below the upper bounds by the Planck result Ω DM h 2 = 0.1199 ± 0.0027 [179] in combination with the WMAP data [180] (with a 10% theoretical uncertainty).Besides, we require that the Spin-Independent (SI) and Spin-Dependent (SD) DM direct detection constraints, for example, the LUX [181], XENON1T [182,183], and PandaX-4T [184,185]-should be satisfied.We should note that such DM constraints can in fact be relaxed in gauge mediation scenarios, as the gravitino mass can be as light as O(1)eV, which can account for the present DM relic abundance for m 3/2 ∼ keV.In our conservative survey, such DM constraints are included in our numerical results.
(v) The muon g − 2 anomaly needs to be explained within the 2σ range for the recent updated combined data.We require the new physics contributions to ∆a µ to lie within We use NMSSMTools_5.6.2 [109,186] to numerically scan the parameter spaces for muon the g − 2 explanation with the following low-energy NMSSM inputs Instead of using the soft SUSY breaking parameters m 2 Hu , m 2 H d and m 2 S , one usually trades them for m Z , tan β and µ ef f = λs by implementing the scalar potential minimization conditions.The squarks and gluino are chosen to be heavy to evade the stringent LHC constraints on the colored sparticles.The bias terms that trigger the collapse of DWs are related to the explicit Z 3 breaking χm 3/2 parameter.We check that the presence of small explicitly Z 3 violation bias terms with χm 3/2 ∼ 10 −12 eV will not alter the Z 3 invariant NMSSM spectrum.
The present-day relic abundance of stochastic GW backgrounds can be calculated numerically with the tensions of the DWs and the related bias terms for each parameter point that satisfy the constraints from (i) to (v).The BBN and t dom constraints in (4.11) are also taken into account.
We have the following discussions for our numerical results • The χS 2 case: From our estimation, we adopt the following range of χm 3/2 10 −16 eV < χm 3/2 < 10 −12 eV, (4.33) in our numerical scan.The peak amplitudes of the GW backgrounds at the present time versus the corresponding peak frequencies are shown in the left panel of Fig. 1.The corresponding GW spectra (with the NANOGrav observation data denoted by the red contour) are shown in the right panel of Fig. 1.We can see that the GW spectrum for many of the survived parameter points can have overlap with the reconstructed posterior distributions for the NANOGrav observations of nHz GWs so as that the DW explanation from Z 3 -NMSSM is viable.
We show in the left panel of Fig. 2 our numerical results on the tensions of the DWs versus their collapse time.The parameter points, whose GWs spectra can overlap with the NANOGrav signal region, are marked with the color blue.From the distributions of the blue points, we can see that the PTA data can impose stringent constraints for the muon g − 2 favored parameter space.As discussed in the previous section, the lower bounds from σ (hence giving lower bounds for m sof t ) can exclude some regions that can explain the muon g − 2 with smaller m sof t .From the right panel of Fig. 2, we can see that the range of χm 3/2 that can account for the NANOGrav nHZ GW background data should lie within 10 −15 eV ∼ 10 −13 eV.As anticipated, there is approximate linear dependence between χm 3/2 and f peak .
We show in the left panel of Fig. 3 the SUSY contributions to muon g − 2, the values ∆a µ , versus the present-day peak amplitudes of stochastic GW background Ω gw h 2 for those survived points whose GWs spectra can overlap with the NANOGrav signal region.We can see that there are still large parameter regions that can explain the muon g − 2 anomaly within 1σ range.Large values of Ω gw h 2 always favor large σ, hence large m sof t .As the SUSY explanation of ∆a µ in general favor small m sof t of order O(10 2 )GeV, the preferred ranges of ∆a µ that can account for the muon g − 2 anomaly favor low Ω gw h 2 .The ranges of the λ and κ parameters of NMSSM that can account for the NANOGrav GW background observations are shown in the right panel of Fig. 3.Both of them are upper bounded to be be less than 0.18.
We can have similar discussions for this case.The peak amplitudes of the GW backgrounds at the present time versus the corresponding peak frequencies is shown in the left panel of Fig. 4. The corresponding GW spectra are shown in the right panel of Fig. 4. Again, the GW spectrum for many the surviving parameter points can have overlap with the reconstructed posterior distributions for the NANOGrav observations of nHz GWs.So, it is obvious that the collapse of DWs from approximate   In the left panel of Fig. 5, we show our numerical results on the tensions of the DWs versus their collapse time.From the distributions of the blue points, which denote those parameter points whose GWs spectrum can be consistent with the NANOGrav signal region, we can see that the PTA data can again impose stringent constraints for the muon g − 2 favored parameter space.Similarly, the lower bounds from σ that lead to lower bounds for m sof t can exclude some regions that can explain the muon g − 2 with light sparticle masses.From the right panel of Fig. 5, we can see that the range of χm 3/2 that can account for the NANOGrav nHZ GW backgroun5d data should lie within 10 −11 eV ∼ 10 −9 eV.The correlations between χm 3/2 and f peak are not obviously linear, as the bias terms in this case can differ in a relatively wide range, which involves v 2 ,tan β and m sof t .We show in the left panel of Fig. 6 the total SUSY contributions to muon g − 2, the value ∆a µ , versus the present-day peak amplitudes of stochastic GW background Ω gw h 2 for those survived points whose GWs spectra can overlap with the NANOGrav observations.Again, there are still large parameter regions that can explain the muon g − 2 anomaly within 1σ range.The NANOGrav GW background data can set upper bounds for λ and κ parameters of NMSSM, which should be less than 0.14 (see the right panel of Fig. 6).

Conclusions
The spontaneous breaking of the discrete Z 3 symmetry in approximate Z 3 -invariant Nextto Minimal Supersymmetric Standard Model by the VEVs of Higgs fields can lead to the formation of DWs in the Universe.Such potentially problematic DWs can collapse and lead to the stochastic GW background signals observed by PTA collaborations due to some explicitly Z 3 breaking terms in the NMSSM effective superpotential and scalar potential.
In the presence of a hidden sector, such terms may originate from the geometric superconformal breaking with holomorphic quadratic correction to frame function when the global scale-invariant superpotential is naturally embedded into the canonical superconformal supergravity models.The smallness of such Z 3 breaking mass parameters in the NMSSM may be traced back to the original superconformal invariance.
Naive estimations indicate that SUSY explanation to the muon g − 2 anomaly can have tension with the constraints on SUSY by PTA data, because large SUSY contributions to ∆a µ in general need relatively light superpartners of order O(10 2 )GeV while the lower bounds for present Ω 0 gw can also set the lower bounds for the tension of DWs σ, leading approximately to O(1)TeV lower bound for m sof t with the BBN upper bounds for the decay time of DWs.We calculate numerically the signatures of GWs produced from the collapse of DWs by the corresponding bias terms.We find that the observed nHZ stochastic GW background by NANOGrav, etc., can indeed be explained with proper tiny values of χm 3/2 ∼ 10 −14 eV for χS 2 case (and χm 3/2 ∼ 10 −10 eV for χH u H d case), respectively.Such tiny values can be natural consequences of higher-dimensional operators in the Kahler potential with U (1) R symmetry, which do not have the notorious tadpole problem.Besides, there are still some parameter points-whose GWs spectra intersect with the NANOGrav signal region, can explain the muon g − 2 anomaly to 1σ range.

Figure 1 .
Figure 1.The left panel shows the present-day peak amplitude of stochastic GW background Ω gw h 2 versus the peak frequency f peak for the parameter points of χS 2 case.The corresponding GW spectra (with the NANOGrav observation data denoted by the red contour) are shown in the right panel.Each parameter point satisfies the experimental constraints (i)-(v).

Figure 2 .
Figure 2. The tensions of the DWs versus their collapse time for the parameter points of χS 2 case are shown in the left panel.The parameter points that can account for the NANOGrav GW background observations are marked with blue color.Each parameter point satisfies the experimental constraints (i)-(v).The values of χm 3/2 versus the peak frequencies are shown in the right panel.

Figure 3 .
Figure 3.The tensions of the DWs versus their collapse time for the parameter points of χS 2 case are shown in the left panel.The parameter points that can account for the NANOGrav GW background observations are marked with the color blue.Each parameter point satisfies the experimental constraints (i)-(v).The values of χm 3/2 versus the peak frequencies are shown in the right panel.

Figure 4 .
Figure 4.The left panel shows the present-day peak amplitude of stochastic GW background Ω gw h 2 versus the peak frequency f peak for the parameter points of χH u H d case.The corresponding GW spectra (with the NANOGrav observation data denoted by the red contour) are shown in the right panel.Each parameter point satisfies the experimental constraints (i)-(v).

Z 3 -
invariant NMSSM with the small explicitly Z 3 breaking χH u H d term (and its soft SUSY breaking B-terms) can explain the PTA observations.

Figure 5 .
Figure 5.The tensions of the DWs versus their collapse time for the parameter points of χH u H d case are shown in the left panel.The parameter points that can account for the NANOGrav GW background observations are marked with the color blue.Each parameter point satisfies the experimental constraints (i)-(v).The values of χm 3/2 versus the peak frequencies are shown in the right panel.

Figure 6 .
Figure 6.The tensions of the DWs versus their collapse time for the parameter points of χH u H d case are shown in the left panel.The parameter points that can account for the NANOGrav GW background observations are marked with the color blue.Each parameter point satisfies the experimental constraints (i)-(v).The values of χm 3/2 versus the peak frequencies are shown in the right panel.

Table 1 .
The tension, energy density of DWs from Z 3 -invariant NMSSM and the SUSY contributions to ∆a µ for some benchmark points that can also satisfy all the low-energy collider constraints.