Investigating multiple solutions to boundary value problems in constrained minimal and non-minimal SUSY models

We investigate the physical origins of multiple solutions to boundary value problems in the fully constrained MSSM and NMSSM. We derive mathematical criteria that formulate circumstances under which multiple solutions can appear. Finally, we study the validity of the exclusion of the CMSSM in the presence of multiple solutions.


Introduction
After the discovery of the Higgs boson with a mass of M h = (125.10 ± 0.14) GeV [1][2][3] and the non-discovery of supersymmetric (SUSY) particles at the LHC, it becomes clearer that pure weak-scale supersymmetry may not be realized in nature. Although the general Minimal Supersymmetric Standard Model (MSSM) may be difficult to fully exclude, the constrained MSSM (CMSSM), which is inspired by minimal supergravity, is in tension with precision observables at more than 90% confidence level [4,5]. One reason for the tension is that in the CMSSM, by construction, all sfermion masses are of the same order. In this case, however, observables such as the Higgs boson mass, Dark Matter relic density and the anomalous magnetic moment of the muon cannot be explained simultaneously by the MSSM, because M h ≈ 125 GeV requires multi-TeV stops, while other observables like (g − 2) μ prefer sub-TeV sleptons.
However, in Refs. [6][7][8] it was discovered that there may be multiple MSSM parameter sets which fulfill the same CMSSM boundary conditions. The mathematical reason for this phenomenon is that the CMSSM is formulated as a boundary value problem (BVP), where the running MSSM This research was supported by the German DFG Research Unit New Physics at the LHC (FOR 2239) and the German DFG Excellence Strategy Quantum Universe (EXC 2121). a e-mail: daniel.meuser@desy.de b e-mail: alexander.voigt@physik.rwth-aachen.de DR parameters are fixed by input values at different renormalization scales. The parameters at the different scales are connected via a set of differential equations, the so-called renormalization group equations (RGEs). Formally, such a BVP may have no, one, or multiple solutions for the MSSM parameters.
In order to make a statement about the validity of the CMSSM, all possible solutions to the BVP must be studied. However, the BVP solving algorithm used in the global fitting analyses of Refs. [4,5] can at most find one solution and may miss further ones. This raises the question to which extent the CMSSM remains challenged in the presence of multiple solutions of the BVP.
In the present paper we systematically study the physical origin of the multiple solutions in the CMSSM. In doing this, we go beyond the scope of Refs. [6,7] and derive mathematical criteria that formulate circumstances under which multiple solutions can appear. In addition we study the influence of the chosen low-energy observables that fix the electroweak gauge couplings, which appear to play in important role for the occurrence and the number of multiple solutions. Next, we apply our newly gained insights to the results presented in Ref. [4] and investigate the influence of multiple solutions on the global fit performed therein. Finally, we extend our analysis to the fully constrained Next-to-Minimal Supersymmetric Standard Model (CNMSSM) and demonstrate that multiple solutions can also occur in constrained non-minimal SUSY models.

CMSSM boundary conditions
The CMSSM is formulated as a BVP, where the DR parameters are fixed at three different scales, see Fig. 1. At the electroweak scale Q = M Z , the DR gauge and Yukawa At the SUSY scale Q = M S ≡ √ mt 1 mt 2 , where mt i denotes the ith DR stop mass, the parameters |μ| and Bμ are fixed by the two electroweak symmetry breaking (EWSB) equations. This leaves the following five free parameters of the CMSSM:

CNMSSM boundary conditions
In the Z 3 symmetric NMSSM the parameters μ and Bμ are absent and are replaced by the new parameters λ, κ, v s , m 2 s , A λ and A κ [9]. In the following study we consider the fully constrained Z 3 symmetric NMSSM (CNMSSM) [10][11][12][13][14], where all soft-breaking parameters are unified at the GUT scale as in Eq. (1) with the additional constraints see Fig. 2. This leaves seven free parameters t β , m 2 0 , M 1/2 , A 0 , λ, κ, and v s , of which three are fixed by the NMSSM EWSB equations. Since λ and κ are dimensionless parameters that enter all NMSSM β functions at sufficiently high order, we take them as input at Q = M S and fix the parameters m 2 0 ,

Fig. 2 CNMSSM boundary value problem
M 1/2 and A 0 by the three EWSB equations via the semianalytic approach, see below. As a result, we are left with the following four free CNMSSM parameters: In the following we trade v s for μ eff = λv s / √ 2.

Matching to low-energy observables
For the study of multiple solutions in the C(N)MSSM, the calculation of the DR electroweak gauge couplings g 1 (Q) and g 2 (Q) is of particular importance. In our analysis we use FlexibleSUSY 2.1.0 [15,16], where the g i (Q) are determined from the DR electromagnetic fine-structure constant α(Q) and the DR weak mixing angle θ W (Q) as Since version 2.0.0, FlexibleSUSY offers the following two possibilities to calculate θ W (Q): 1. The Z 0 pole mass M Z and the Fermi constant G F can be chosen as input, in which case the weak mixing angle is calculated as where Δr (Q) is a function of the one-loop Z 0 and the W ± self-energies [17]. 2. The Z 0 pole mass M Z and the W ± pole mass M W can be chosen as input, in which case the weak mixing angle is calculated as In Eq. (7) m W (Q) and m Z (Q) denote the DR W ± and Z 0 masses, respectively, calculated as with V ∈ {W, Z } and Π V being the corresponding oneloop vector boson self-energy.
Both Eqs. (6) and (7) are equivalent at the one-loop level, but differ at higher orders. The difference involves in particular products of the one-loop vector boson self-energies. This difference is of crucial importance for the existence of multiple solutions at small values of |μ|, as is shown in Sect. 3.1.

Boundary value problem solvers
The most common way to solve the CMSSM BVP is the "two-scale solver" (TSS), also known as "running and matching". In this approach the spectrum generator numerically integrates the RGEs and imposes the boundary conditions at each scale in an iterative way. If the iteration converges, the algorithm has found one solution of the BVP. In the CMSSM one can use the TSS to choose m 2 0 , M 1/2 , A 0 , and sign(μ) as input and the parameters μ and Bμ become an output at the SUSY scale M S . As an illustration we show in Fig. 3 the relation between m 2 0 and μ for the parameter set [7] t β (M Z ) = 40, M 1/2 = 660 GeV, A 0 = 0, and for both signs of μ. The blue dashed line denotes all points where the TSS could find a solution. We see for instance that for m 0 = 2 TeV the TSS could find a solution for sign(μ) = +1, but not for sign(μ) = −1. For m 0 = 3 TeV the TSS finds one solution for each sign(μ) and for m 0 3.5 TeV no solution is found. The CNMSSM BVP, however, cannot be solved in this way with the TSS. There, if one chooses m 2 0 , M 1/2 , and A 0 as input at the GUT scale, the parameters λ, κ and v s would need to be fixed by the EWSB equations at the scale M S . Since λ and κ are dimensionless parameters which enter most NMSSM β functions, the iteration between the scales becomes unstable and the TSS does not converge.
This downside of the TSS led to the development of the semi-analytic solver (SAS) [16,[18][19][20][21][22], which allows one to exchange the role of input and output parameters in constrained models. To achieve that, the SAS exploits the general structure of the RGEs to decompose the soft-breaking parameters in terms of GUT parameters and dimensionless coefficients c The coefficients c ( j) i are scale dependent but independent of the (dimensionful) GUT parameters. They are determined numerically by solving the BVP for different values of m 2 0 , M 1/2 , A 0 , and Bμ(M X ) with the TSS. Once the coefficients are known, Eq. (10) can be solved for the soft-breaking parameters, which become an output of the algorithm.
The property of the SAS to exchange input and output parameters has several advantages over the TSS: -If the relation between input and output parameters is not injective (which occurs in both the CMSSM and CNMSSM), scanning over one parameter while obtaining the other one as output, and vice versa, allows the search for multiple solutions of the BVP. This procedure is shown by the red solid line in Fig. 3, where μ(M S ) is used as input and m 2 0 is output. One immediately sees that with the SAS one can find up to four solutions for μ around m 0 ≈ 3.3 TeV, which have not been found by the TSS. In Sect. 3 we use this procedure to study in depth the physical origin of the multiple solutions found in Refs. [6][7][8].
-In models where the TSS would require dimensionless parameters to be output, as for example in the CNMSSM or CE 6 SSM [19,23], the SAS enables one to take them as input, which yields a stable iteration between the high and low scales. In Sect. 4 we use this feature in the CNMSSM to take the parameters λ(M S ), κ(M S ), and v s (M S ) as input and obtain m 2 0 , M 1/2 , and A 0 as output.
On the other hand, there are scenarios where it is advantageous to use the TSS: -If the derivative μ (m 0 ) approaches zero, as happens for example in Fig. 3 around μ ≈ 800 GeV, scanning over μ is no longer suitable. Invoking the TSS to vary m 0 instead and receiving μ as an output allows to study a region of parameter space which cannot be accessed via the SAS as easily. -In cases where both solvers find the same unique solution it is a priori not clear whether the TSS or the SAS is the better choice. A combination of both allows a more complete study and a comparison validates the equivalence of the two solvers. In those cases, the TSS in general converges much faster than the SAS does.
In the following section we use the semi-analytic approach to systematically search for multiple solutions to the BVP of the CMSSM and study their origin in depth. For our analysis we use the SAS implemented in FlexibleSUSY 2.1.0. [15,16].

Effects from light SUSY particles
In this section we study the occurrence of multiple solutions in the CMSSM for small values of |μ| (see Fig. 3), which were first observed in Refs. [6,7]. Without limitation of generality we restrict our analysis to positive values of μ by making use of the approximate mirror symmetry m 0 (μ) ≈ m 0 (−μ), which allows two values for μ for fixed m 2 0 as long as |μ| is not too large. As was shown in Refs. [6,7], even for fixed sign(μ) the function m 0 (μ) is not necessarily bijective: In the region 0 < |μ| ≤ 100 GeV the function has several turning points where the derivative m 0 (μ) exhibits singular behavior, see The origin of these kinks can be traced back to singular chargino and neutralino contributions to the one-loop vector boson self-energies Π Z and Π W (see Fig. 6), which are used to calculate the electroweak gauge couplings at the scale Q = M Z , as described in Sect. 2.1. Since the gauge couplings contribute to every β function of the MSSM, a singular point in Π V (μ) translates into a singular point in m 0 (μ). Note that the precise dependence of the g i (Q) on the Π V depends on The vector boson self-energies Π V depend on μ via the chargino and neutralino masses, which enter the one-loop diagrams shown in Fig. 6. Expressed in terms the loop functions H 0 and B 0 [17], each diagram gives a contribution to the self-energies. In Eq. (11) p 2 denotes the external momentum squared, m 1,2 are the running masses of the fermions in the loop, and C L and C R are vertex coefficients derived from the Lagrangian. All running quantities are evaluated at the renormalization scale Q = M Z . In total there are fourteen electroweakino diagrams contributing to Π Z and eight to Π W . Of these, four (Z 0 ) and two (W ± ) contain only higgsino-like charginos and neutralinos, whose masses are approximately given by |μ|. Therefore, only those six diagrams are responsible for the kinks. The singularities in Π V (μ) now have the following deeper origin: In the DR renormalization scheme the finite part of the B 0 function appearing in Eq. (11) is given by where In Fig. 7 the real and imaginary part of B 0 are shown exemplary as a function of the mass m 2 for fixed m 1 . One finds that as soon as m 2 becomes small enough such that both particles in the loop go on-shell simultaneously, B 0 acquires an imaginary part and the derivatives of the real and the imaginary part are singular for this value of m 2 . This is due to the square root function in Eq. (13b) not being differentiable when its argument vanishes. The radicand, which is just the Källén function 2 λ( p 2 , m 2 1 , m 2 2 ), has two roots of which the one with positive sign is located at m 2 ≈ 41 GeV and the one with negative sign is located at m 2 ≈ 141 GeV in Fig. 7. Both real and imaginary part are not differentiable at the zero m 2 ≈ 41 GeV, but they are at the other due to a cancellation between the f B -terms in Eq. (12). The appearance of an imaginary part is in accordance with the optical With this knowledge one would expect four higgsino-like (and thus μ-dependent) combinations (χ 0 to give rise to singularities in the Z 0 self-energy, and two combinations (χ 0 1 χ ± 1 , χ 0 2 χ ± 1 ) for the W ± self-energy. Following, there should be four spikes for the input parameters {M Z , G F } and six for {M Z , M W } in total. However, there are only two and four spikes in m 0 (μ), respectively. Evaluating Eq. (11) in the limit of a vanishing Källén function renders it in a form in which the disappearance of some singularities becomes manifest. To this end we use the identities where we have neglected constants as well as terms which are proportional to the one-point function A 0 [17]. Plugging Eq. (15) in (11) and setting λ = 0, we arrive at We plug the root Eq. (14) with positive sign into Eq. (16) to obtain Equation (17) is the non-vanishing part of a self-energy graph in the limit where the Källén function λ( p 2 , m 2 1 , m 2 2 ) is zero. From this we infer that whenever C L + C R is zero in the limit λ → 0, the corresponding diagram does not contribute a singularity to Π V .
Since the entries of the chargino mixing matrices U − and U + are independent, there is no relation which would guarantee the cancellation C L + C R = 0, as long as at least one chargino appears in the loop. Hence, the diagrams with χ 0 1 χ ± 1 , χ 0 2 χ ± 1 , and χ ± 1 χ ± 1 in the loop contribute nonvanishing singularities to Π V .
The other three relevant diagrams are pure χ 0 i χ 0 j contributions to the Z 0 self-energy. For a Z 0 χ 0 i χ 0 j vertex the relation C * L = −C R holds [17]. Since the C L/R are in general complex quantities, this property is not sufficient for C L + C R to vanish. If the neutralinos are identical, however, C L and C R become real and cancel when added.
We conclude that the Z 0 self-energy diagrams with identical neutralinos in the loop do never give a spike, since the singular terms in Π Z have a vanishing coefficient in the limit λ → 0. The singularities (kinks or spikes) originate only from diagrams with light χ 0 Note, that at the singular points spikes or kinks can appear in the overall vector boson self-energies. However, only spikes lead to multiple solutions, because there the sign of the derivative around the singular point changes, which results in a turning point. Whether a singularity causes a spike or a kink depends on the relative sign between the B 0 function which causes the singularity and other loop corrections to the self-energy. The relative signs generally depend on the regarded model as well as on the input parameters.

Effects from non-linear parameter inter-dependencies
In Refs. [6,7] another CMSSM parameter region with multiple solutions was found around μ ≈ −500 GeV, see Fig. 8. In this section we investigate this parameter region with multiple solutions, which we find are not caused by light particles in vector boson self-energies, but by intricate nonlinear parameter inter-dependencies.
The non-monotonic behavior of m 0 (μ) in the region μ ≈ −500 GeV is depicted in Fig. 8. Multiple solutions appear in this region, because the function m 0 (μ) changes its sign such that m 0 (μ) has a minimum at μ = −540.5 GeV. For too small values of μ −545 GeV there is no physical solution because the running masses of the neutral Higgs bosons become tachyonic in this region of parameter space. In the following we study the parameter interplay which is responsible for the existence of the minimum.
To this end, we derive an approximate relation between m 0 and μ. Our starting point is the tree-level EWSB equation which is imposed at the SUSY scale M S . Hence, all renormalization scale dependent quantities are evaluated at M S . For our set of input parameters, t 2 β 1 holds and we can also neglect the m Z -term, which implies To relate those parameters to m 2 0 , we make use of the RGEs At leading order the β functions are independent of μ and fulfill the hierarchy β m 2 hu > β m 2 h d . This hierarchy arises because the up-type β function is proportional to the squared top-Yukawa coupling whereas the down-type one contains only down-type fermion contributions. It causes m 2 h d to be close to m 2 0 and at the same time allows m 2 h u to become negative -which is usually necessary for EWSB to occur. We will therefore only replace m 2 h d in Eq. (19) by (20) and arrive at This, to leading order, implies that a minimum of m 2 0 appears when For our choice of parameters, a numerical analysis leads to so Eq. (22) is fulfilled for μ = −540.5 GeV. A comparison with Fig. 8 validates that this indeed is the position of the minimum around which multiple solutions occur. For differently chosen input parameters the existence of a point where the derivative of m 2 h u is of the appropriate size to create a minimum is not ensured. An example for this is given in Sect. 3.4.

Multiple solutions for M 1/2 and A 0
The SAS allows one to exchange input and output parameters of constrained (SUSY) models. This property was used in the previous sections to exchange the role of m 2 0 and μ in the CMSSM to examine the function m 2 0 (μ) for (non-)injectivity. However, the SAS is not restricted to the exchange m 2 0 ↔ μ (SAS1) and further choices are possible, see Table 1. In the following we apply the SAS to study the relations M 1/2 ↔ μ (SAS2) and A 0 ↔ μ (SAS3) to investigate the implications of their non-bijectivity.
The semi-analytic solver can treat the GUT parameter M 1/2 as output and μ as input (SAS2) by solving the EWSB equations for M 1/2 . To do so, one makes use of the semianalytic ansatz Eq. (10) again. The equation involving the m h i is quadratic in M 1/2 , in contrast to being linear in m 2 0 , which allows for up to two distinct solutions for fixed μ. In Fig. 9a, b we scan over μ (SAS2) and M 1/2 (TSS) for A 0 = 0, t β = 40 and two different values m 0 ∈ {3000, 3500} GeV. As expected, the semi-analytic solver finds two branches, one for each solution of the quadratic equation. The TSS only finds solutions which are obtained from SAS2 as well. The solutions at M 1/2 = 660 GeV (vertical black lines) are the same as found in Sect. 3.1 using SAS1. Thus, scanning over μ while using M 1/2 as an output parameter does not give new solutions compared to the case where m 2 0 is output. Similarly to SAS2, A 0 can be treated as output parameter by solving the EWSB equations for A 0 (SAS3) instead. Also in this case the ansatz is quadratic in A 0 , so up to two distinct solutions are possible. In Fig. 10a Fig. 3 A 0 (μ) and μ(A 0 ) for M 1/2 = 660 GeV, t β = 40, and m 0 ∈ {3000, 3500} GeV. We find that for m 0 = 3500 GeV (Fig. 10b) the SAS can scan over both solution branches, while for m 0 = 3000 GeV (Fig. 10a) the two branches have merged. The multiple solutions around μ ≈ 0 have vanished and a region without solutions (|μ| 300 GeV) has emerged. The curve μ(A 0 ) becomes flat around A 0 ≈ 1000 GeV and SAS3 is not suitable to find solutions in the region where the two solution branches merge; the two-scale solver, however, is able to find more solutions in this case. Besides this, the multiple solutions that we obtain for SAS3 are again the same as found for SAS1 in Sect. 3.1.
In conclusion, scanning over μ while receiving either M 1/2 or A 0 as output does not give us any new solutions compared to the m 2 0 output search strategy. We nevertheless have been able to reproduce the previous results in a consistent manner and also our expectation of finding two branches of solutions for parameters of mass dimension one has been   Fig. 3 fulfilled. Furthermore, the same singular structures for small values of μ have appeared for SAS2 and SAS3. Note also, that both the two-scale and the semi-analytic solver had to be used to find all solutions in the considered parameter regions.

Is the CMSSM still challenged?
In this section we investigate the relevance of multiple solutions for globally fitting the CMSSM to experimental and observational data. As a reference point we use Ref. [4] ("Killing the CMSSM softly"), in which the program Fittino was used. The idea was to scan over a reasonable region of the CMSSM input parameter space and to determine the fit point which is the most compatible with the considered observables. Furthermore, this paper was the first to derive a consistent p-value for the CMSSM from toy experiment.

The results of "Killing the CMSSM Softly"
We briefly list the observables which have been used in the analysis of Ref. [4]. As for the precision observables there are the anomalous magnetic moment of the muon a μ , the effective weak mixing-angle sin θ eff , the top quark and W boson masses as well as the b quark/B meson branching ratios. Additionally, different combinations of Higgs observables, like e.g. the SM Higgs boson mass and its decay channels, as observed by the ATLAS/CMS experiments at the LHC, and the dark matter relic density Ωh 2 , as measured by the Planck collaboration, are incorporated.
The following parameter values were found to give the best accordance between measured and predicted observables and are hereinafter referred to as "best-fit point" parameters: It should be noted that these values were determined with the spectrum generators SPheno 3.2.4 [24,25] and FeynHiggs 2.10.1 [26][27][28][29][30][31][32][33][34][35][36], whereas the following analysis is based on the predictions of The chargino masses are determined by M 2 and μ. The lightest chargino is thus able to escape the LEP bound m χ ± 1 > 94 GeV [37]. The lightest supersymmetric particle (LSP) is the lightest neutralino with mass ∼ 400 GeV and provides a dark matter candidate. At the best-fit point the predicted dark matter relic density Ωh 2 is in agreement with the experimental observations. The Standard Model prediction for the anomalous magnetic moment of the muon a μ deviates from the observed value at a 3.5σ level. To account for this, a successful SUSY model is expected to give an additional contribution to a μ of the order 30 × 10 −10 [38]. The best-fit point predicts a correction a SUSY μ ∼ 4 × 10 −10 , which is far too small. This discrepancy is the main reason for the smallness of the p value; when taking (g − 2) μ into account, one finds p = (4.9 ± 0.7)%, whereas leaving it out results in the much larger p = (51 ± 3)%.
When performing a global fit of a model to known observables, all possible mathematical solutions to the formulated boundary value problem need to be taken into account. The TSS of SPheno, however, which was used in Ref. [4], does not necessarily find all solutions. For this reason we study in the following section whether further solutions to the CMSSM BVP can be found around the best-fit point with the semi-analytic approach.

Multiple solutions around the best-fit point
In this section we perform an analysis similar to the one of Sect. 3.3 for the best-fit point [cf. Eq. (24)]. All except one of the free CMSSM GUT input parameters are set to their bestfit values and the remaining parameter is varied subsequently with the TSS of FlexibleSUSY. To find possible multiple solutions we repeat the same procedure but scan over μ(M S ) by means of the semi-analytic solvers SAS1-SAS3. The results are shown in Fig. 11.
As for Fig. 11a we recognize the same overall relation between m 0 and μ as in Fig. 3 and the multiple solutions around μ = 0 appear as well. At the lower end of the curve, the multiple solutions due to non-linear parameter inter-dependencies are non-existent. There, for instance, and so Eq. (22) is not fulfilled. In Fig. 11b, c the semi-analytic solver does not find any additional solutions compared to the two-scale solver. In the case of Fig. 11c the curve μ(A 0 ) becomes too flat around A 0 ≈ 1500 GeV to allow an efficient scan over μ with the SAS3. For some values of the scan parameters, for instance the region where 0 < M 1/2 < 400 GeV in Fig. 11b, we do not find a physical solution due to either tachyonic down-type sleptons, up-type squarks, or Higgs bosons.
Concluding, we find that the best-fit point is far off from the regions in which multiple solutions can occur. One reason is that around the best-fit point all SUSY particles are heavy enough to escape the experimental constraints, while our previous analysis has shown that additional solutions tend to occur in regions of parameter space where at least some superpartners become light.

Multiple solutions in the CNMSSM
In this section we study the occurrence of multiple solutions for a given set of parameters {m 2 0 , λ, κ} in the CNMSSM. Exchanging m 2 0 with μ eff allows us to search for multiple solutions in a similar fashion as we did for the CMSSM. It has to be noted, however, that we now also take the GUT parameters M 1/2 and A 0 to be output of our algorithm and not input as we did in the CMSSM. The CMSSM GUT parameters and μ are varied around the best-fit point from Ref. [4] using the two-scale solver and the different semi-analytic solvers (cf. Table 1). In a the best-fit point lies in a region without multiple solutions. In b, c no multiple solutions appear at allapart from the expected ones which differ by sign(μ) 4.1 Study of multiple solutions with the semi-analytic approach As described in Sect. 2.2, the CNMSSM is formulated as a BVP with the universal GUT parameters {m 2 0 , M 1/2 , A 0 }. In order to study multiple solutions in the CNMSSM, we make use of the semi-analytic equations, which allows us to take the SUSY scale parameters {λ, κ, μ eff } as input. Specifying the dimensionless quantities λ and κ yields a good stability of the underlying solving algorithm. μ eff can now be used as a scan parameter as was done for the CMSSM. As a result, we formulate the CNMSSM BVP in terms of the input parameters and obtain {m 2 0 , M 1/2 , A 0 } as output. In Fig. 12 we show a scan over μ eff (M S ) for fixed values of t β , λ, and κ. The output parameters {m 2 0 , M 1/2 , A 0 } are shown on the abscissae. The curves are discontinuous due to the existence of two distinct solutions branches, which differ from each other in their sign of M 1/2 . In order to distinguish the branches, solutions with sign(M 1/2 ) = ±1 are marked as red and cyan dots, respectively. If we were able to pre-select one of the branches, a scan should yield a continuous relation between the output parameters and μ eff .
In contrast to the CMSSM, where we found multiple solutions, now we do not find any physical solutions in the parameter region |μ eff | M Z ,W at all. The reason for this is that, in our scenario, as μ eff tends to 0, the dimensionful output parameters become small and tachyons appear in the particle spectrum, i.e. the solutions of the BVP have to be discarded. Furthermore we find a highly non-linear dependence of m 2 0 on μ eff for each sign(M 1/2 ), see Fig. 12a, which results in up to three solutions around m 2 0 0. This non-linear dependence can be understood as follows: When determining the three GUT parameters from the three EWSB equations of the CNMSSM, the strongest restriction to the value of m 2 0 comes from the equation [9] For small λ and κ, the parameters m 2 s and A κ are approximately constant between the GUT and the EW scale [ From Fig. 12c we can also see A 0 = a(v s )v s with some positive function a(v s ), which depends only weakly on v s , and so for κ < 0. The term with positive sign is only approximately quadratic in v s and takes the shape of a slightly tilted parabola. It is the sum of both terms which causes the function m 2 0 (μ eff ) to behave in the observed non-monotonic way. The relation between M 1/2 and μ eff is nearly linear, as one would expect for two parameters of mass dimension one. The two branches with different sign(M 1/2 ) are related to one another by an approximate central symmetry, see Fig. 12b. The reason for the resulting discontinuous cross-like shape is that, as explained above, for a given value of μ eff the semianalytic solver usually can find a solution for one value of sign(M 1/2 ), but not simultaneously for the opposite one. A similar behavior can be found for A 0 (μ eff ) in Fig. 12c, where the semi-analytic solver can find one solution for fixed μ eff at most. Here, however, the sign of A 0 is fully determined by the sign of μ eff and the output parameter differs only slightly between the two solution branches.

Mass spectra for two different solutions of a single CNMSSM parameter point
In order to see the importance of multiple solutions in the CNMSSM, we compare the mass spectra of two different solutions for the parameter point The two solutions share the same m 2 0 , but have different M 1/2 , A 0 and μ eff , see Table 2 and Fig. 13. As a result, the points have different pole mass spectra as shown Fig. 13b. By comparing the predicted Higgs boson pole masses with the experimentally measured value of M h = (125.10 ± 0.14) GeV, point 2 can be excluded, while point 1 may still be taken into consideration. However, none of the points can correctly predict the observed Dark Matter relic density because the lightest stau is the lightest supersymmetric particle (LSP). This is in agreement with the analysis of Ref. [14], finding that the condition A 0 ∼ −M 1/2 /4 must be fulfilled in order for a CNMSSM parameter point to produce the observed dark matter relic density.
In any case, our analysis shows that the different solutions can have a significantly different phenomenology and the SAS is a useful tool to find and study multiple solutions in this model. However, since M 1/2 and A 0 are output parameters in our SAS formulation of the CNMSSM BVP, viability conditions such as the ones given in Ref. [14] cannot (a) (b) Fig. 13 Two parameter points fulfilling the CNMSSM boundary condition given in Eq. (31). The points differ in their values for μ eff , M 1/2 , and A 0 (see Table 2), which results in two vastly distinct pole mass spectra shown in (b) be enforced from the start and have to be tested on the output parameters.

Conclusions
In Refs. [6][7][8] the appearance of multiple solutions to the BVP of the CMSSM has been discovered and studied. In the present paper we have investigated the deeper origin of these multiple solutions. The study was made possible by the semianalytic BVP solver implemented in FlexibleSUSY, which allows to exchange input and output parameters to search for turning points of inverse functions of the BVP and thus allows the systematic search for multiple solutions. We could trace their appearance back to two phenomena: -Light neutralinos and charginos can lead to singular points in the one-loop W ± and Z 0 self-energies, which translate to singular points in the function m 2 0 (μ). At these points the derivative of m 2 0 (μ) can change its sign, which leads to multiple branches in the inverse functions μ(m 2 0 ). The position of the singular points is given by the light neutralino/chargino masses. The number of singular points depends on their couplings to the W ± and Z 0 bosons and on the formulation used to determine the DR weak mixing angle from physical observables. -A non-linear inter-dependence between the parameters m 2 0 and μ can lead to a minimum of the function m 2 0 (μ), resulting in the appearance of multiple branches in the inverse function μ(m 2 0 ) around that minimum.
Furthermore we found that around the CMSSM best-fit point the solution to the BVP is unique and thus the p value for the CMSSM remains unaltered. It has to be noted, however, that the ongoing "Muon g−2 Experiment" at Fermilab [39] might measure a μ to agree with the SM predictions and render the CMSSM fit perfectly good again by lifting the p value to (51 ± 3)%. Finally we have investigated the appearance of multiple solutions in the CNMSSM. We find that: -Multiple solutions around μ eff M W,Z tend to not occur, because in the limit μ eff → 0 also m 2 0 and other supersymmetry-breaking parameters vanish or become of the order of the electroweak scale, which leads to light or tachyonic scalar particles, i.e. unphysical solutions of the BVP. -For small m 2 0 up to three solutions for μ eff can occur due to a non-linear inter-dependence between these two parameters, imposed by the EWSB equations and the β functions. The different solutions may have significantly different physical spectra because of different values for M 1/2 , A 0 and μ eff .
Concluding, we would like to emphasize that in order to investigate the validity of a constrained SUSY model all possible solutions to the BVP must be studied. In combination with the conventional "running and matching" procedure (TSS), the semi-analytic approach (SAS) is a useful tool to search for additional solutions.