Investigating Multiple Solutions in the Constrained Minimal Supersymmetric Standard Model

Recent work has shown that the Constrained Minimal Supersymmetric Standard Model (CMSSM) can possess several distinct solutions for certain values of its parameters. The extra solutions were not previously found by public supersymmetric spectrum generators because fixed point iteration (the algorithm used by the generators) is unstable in the neighbourhood of these solutions. The existence of the additional solutions calls into question the robustness of exclusion limits derived from collider experiments and cosmological observations upon the CMSSM, because limits were only placed on one of the solutions. Here, we map the CMSSM by exploring its multi-dimensional parameter space using the shooting method, which is not subject to the stability issues which can plague fixed point iteration. We are able to find multiple solutions where in all previous literature only one was found. The multiple solutions are of two distinct classes. One class, close to the border of bad electroweak symmetry breaking, is disfavoured by LEP2 searches for neutralinos and charginos. The other class has sparticles that are heavy enough to evade the LEP2 bounds. Chargino masses may differ by up to around 10% between the different solutions, whereas other sparticle masses differ at the sub-percent level. The prediction for the dark matter relic density can vary by a hundred percent or more between the different solutions, so analyses employing the dark matter constraint are incomplete without their inclusion.


Introduction
It has been recently discovered that for a given set of input parameters for the Constrained Minimal Supersymmetric Standard Model (CMSSM) [1][2][3][4][5], there can exist several distinct solutions for the particle spectrum [6]. This is because the renormalisation group equations (RGEs) connecting the high-and low-scale observables form a boundary value problem, which may admit multiple solutions. A priori, this fact has important consequences for experimental collider searches for new physics. Previous analyses that placed bounds on specific regions of the CMSSM parameter space used spectrum generators which missed the additional solutions. The published limits are therefore optimistic, as they ignore these additional, potentially viable solutions. We present here a study which shows that this initial cause of concern is not necessary for many sparticle searches: the phenomenology of the additional solutions are constrained by current and future experiments. However, there are particular cases (notably, those involving charginos), where care is needed and the additional solutions could provide a loop-hole to exclusion.
It was realised some years ago that the minimal supersymmetric standard model has an LHC inverse problem [7]. The problem is essentially that LHC observables possess degeneracies such that the Lagrangian of the MSSM cannot be uniquely determined from them. It was thought though, that in strongly constrained set-ups like the CMSSM with only a few extra parameters and many potential LHC observables, that there would be no inverse problem (although even with fairly precise edge measurements, one could not always discriminate between minimal gauge mediated supersymmetry breaking patterns and the CMSSM [8]). Multiple solutions though, demonstrate that even the CMSSM has its ambiguities and suffers from the opposite of the LHC inverse problem. While the LHC inverse problem has several different parameter points leading to the same set of LHC observables, multiple solutions have the same CMSSM parameters leading to potentially different LHC observables. One must be careful to be aware of the multiple solutions so as to not rule out the model just because one is ignorant of an additional solution. The same caveat potentially applies to other models of SUSY breaking. This does not reintroduce the inverse problem, as sufficiently accurate observables will discriminate between the different solutions, since they are physically distinct.
The CMSSM is a model of softly broken supersymmetry (SUSY) that imposes many relations amongst the possible terms in the broken MSSM. In particular, at the Grand Unified Theory (GUT) scale, the scalar supersymmetric particles have the same mass m 0 , the gauge fermion supersymmetric particles have the mass M 1/2 and the trilinear scalar couplings (Higgs-sfermion-sfermion) are given by a new parameter A 0 multiplied by the corresponding Standard Model Yukawa matrices. Additionally, the parameter tan β specifies the ratio of the two Higgs' vacuum expectation values at the SUSY breaking scale. The final input to the CMSSM is the sign of the Higgsino mass term, sign(µ). The value of µ is set by requiring the calculated Z 0 mass is equal to the measured value (c.f. Eq. (A.1)). The SUSY particle spectrum at any given scale is then determined by solving the RGEs with boundary conditions at the three scales: GUT, SUSY breaking, and electroweak.
Specifying constraints at more than one scale creates a boundary value problem which, in general and in practice [6], can admit multiple solutions. This is different to the case where a single boundary condition has multiple solutions itself (see, for example the case of mSUGRA in Ref. [9], which gave up to three solutions for tan β). The existence of multiple solutions to the CMSSM RGEs for a given sign of µ have evaded the standard SUSY spectrum generators such as SOFTSUSY [10], ISASUSY [11], SPheno [12], and SUSPECT [13]. The standard procedure for calculating supersymmetric spectra by all the current generators is fixed point iteration. This algorithm is blind to the extra solutions because they correspond to unstable fixed points -the standard methods move away from them by necessity. This fixed point iteration method can be modified, and partially improved, by choosing to scan over one of the dimensions of the RGE and also relax one constraint. For instance, one can scan the value of µ, relax the condition on M Z , and then investigate the locations in the one dimensional scan that give the correct prediction for M Z (exp).
The more general method for finding all possible multiple solutions (at least within the available numerical precision and computing time) is to express the RGEs and the input data as an initial value problem. The Cauchy-Lipschitz theorem [14] states that this initial value problem has a unique solution provided the RGE trajectories are Lipschitz continuous (they are unless poles are encountered). For the CMSSM in the dominant third family approximation 1 , one has 11 high-scale parameters to fix, and then runs the RGEs down to the lower scales to predict 11 quantities, most of which are Standard Model observables. The errors in the low-scale predictions are then used to formulate a root-finding problem: find the high-scale parameters that give negligible error for all predictions. We shall show in this paper how to implement this method using both the multidimensional Newton-Raphson and Broyden root finders. We then apply this general solver to well-studied slices of the CMSSM parameter space, to identify regions with multiple solutions and systematically study their phenomenology.
The paper is organised as follows. In Sec. 2 we review the calculation of the spectrum of the CMSSM, and how it is implemented in publicly available computer programs. Sec. 3 discusses the instability of fixed point iteration, and describes an alternative, the shooting method, which can be used to explore the full solution space of a model defined at multiple scales. The details of our numerics are given in Sec. 4. In Sec. 5 we present a map of multiple solutions in the CMSSM (for a few representative points in parameter space) and discuss the phenomenology of these solutions, giving limits based on collider and astrophysical data. Finally, Sec 6 summarises the work with a discussion on the impact of the multiple solutions on experimental exclusions.

Spectrum calculation in the CMSSM
To compute the observable spectrum of SUSY particles in the CMSSM one must determine the MSSM parameters in the model by solving their RGEs with certain boundary conditions. The high-scale boundary conditions of the CMSSM specify a theoretically motivated input on the soft supersymmetry-breaking mass terms in the Lagrangian, and the low-scale boundary conditions match parameters to known data on Standard Model fermion masses and mixings, gauge boson masses and gauge couplings. Over one hundred coupled, homogeneous, non-linear ordinary differential equations evolve the couplings between the different scales. At the high scale there are a certain number of fixed parameters (due to the GUT unification conditions) and the rest are initially free but must be solved for, such that all low scale constraints are satisfied. In this section we detail these boundary conditions and describe the existing method used to numerically compute solutions.
Boundary conditions are specified at three scales: the low scales M Z and M SUSY , and the high scale M GUT . The latter two scales are part of the set of the free parameters and are themselves determined by the boundary conditions. The conditions at these scales are specified below, in terms of: the Yukawa couplings of the top, bottom, and tau, viz. h t , h b , and h τ , respectively; the 3 MSSM gauge couplings, g i , i ∈ {1, 2, 3}; the various SUSY breaking scalar masses, m ϕ ; the gaugino masses, M i , i ∈ {1, 2, 3}; and the SUSY breaking trilinear scalar couplings for top, bottom, and tau, A t , A b and A τ . Also included are constraints on the parameter µ appearing in the superpotential 2 W ⊃ µĤ 1Ĥ2 , and constraints on m 2 3 , the parameter that mixes the two Higgs doublets in the potential V ⊃ m 2 3 H 2 H 1 . These latter two conditions (Eqs. (2.6),(2.7) below) come from the minimisation of the Higgs potential with respect to the neutral components of H 1 and H 2 . In terms of these variables the boundary conditions are The running parameters in Eqs. (2.1)-(2.12) are in the modified dimensional reduction (DRED) scheme [15]. The 'exp' denotes that the value derives from experimental data. We have labelled the input parameter tan β as tan β(input) to discriminate it from the predicted value run down from the GUT scale. The parameters m b,t,τ (M Z ) and g 1,2,3 (M Z ) are obtained from experimental data, subtracting loops due to sparticles and Standard Model particles. The Standard Model electroweak gauge couplings g 1 (exp) and g 2 (exp) are fixed by the Fermi constant G F and the fine structure constant α. The values of g i (exp) are corrected by one-loop corrections involving sparticles. The parameter v(M Z ) ≈ 246 GeV where v 1 and v 2 are the vacuum expectation values of the neutral components of the Higgs doublets H 1 and H 2 , respectively. The modified DRED Z 0 boson mass squared is fixed by are fixed by the soft SUSY breaking mass parameters for the Higgs fields m 2 H i , i ∈ {1, 2}, as well as by the tadpole contributions t i coming from loops. These tadpole contributions have terms linear in µ(M SUSY ) as well as terms that are logarithmic in it, and so Eq. (2.6) is not a simple quadratic equation for µ(M SUSY ). The symbol Π T ZZ (M Z ) denotes the MSSM self-energy correction to the Z 0 boson mass which can be found in Ref. [16] and I 3 is a 3×3 matrix in family space. For further details, see the SOFTSUSY manual [10].
Spectrum calculators for the MSSM that are currently in the public domain, namely ISASUSY [11], SOFTSUSY [10], SPheno [12], SUSEFLAV [17] and SUSPECT [13], use fixed point iteration to find a self-consistent solution to Eqs. (2.1)-(2.12) (or equations very similar to them). The particular algorithm used by SOFTSUSY, for example, is shown in Fig. 1 This iterative algorithm is an example of fixed point iteration, where one attempts to find a set of values that remain invariant under application of a certain function. In the case at hand the values are the parameters of the MSSM and the function (actually a set of functions, one for each parameter) is one iteration of the loop in Fig. 1. If a solution to the system of MSSM RGEs is found then it is considered the one and only solution. Such an assumption is naïve, and it was shown in Ref. [6] that, due to the boundary value nature of the system, multiple solutions can, and do in fact, exist, depending upon CMSSM parameters. Amongst a set of multiple solutions, many of the M GU T scale parameters (the ones that are set by the lower scale boundary conditions) are different, and so are many low-scale observables. Some are more different than others, however. The Higgs potential parameters µ and m 2 3 can be quite different between the different solutions. µ changes the stop, sbottom and stau mixing. This then changes the top, bottom and tau Yukawa couplings at leading log order (i.e. first order in 1/(16π 2 ) log(M GUT /M Z )), because the measurement of the top, bottom or tau mass would include stop, sbottom or stau radiative corrections, respectively 3 . The fact that these couplings change then changes the RGE trajectories of most of the soft masses. Strictly speaking, these are only thus changed beyond leading log order (since it's a loop effect induced by an effect that is already leading log), but the large separation of scales between the weak and supersymmetry breaking scale can enhance the higher order logs.
Whilst fixed point iteration can in some cases use different starting conditions to find multiple solutions, we shall demonstrate in the next section that there are certain solutions which the algorithm can never converge on. We shall also show that such solutions are sometimes important in the context of CMSSM phenomenology.

RGE Solution Methods and Stability
In this section we first describe the potential non-robustness of fixed point iteration due to instabilities, and then detail a method of solving the multi-boundary problem that does not suffer from stability issues.

Fixed Point Iteration
As mentioned in Sec. 2, current publicly available SUSY spectrum generators use fixed point iteration (FPI) to solve the multi-boundary problem relevant to the CMSSM. FPI operates by making a guess for the unknown high-scale parameters and repeatedly runs the RGEs down and up, applying the constraints until the solution no longer changes above some prescribed level. It has the distinct advantage of being very quick, as often only a few iterations are required in order to achieve convergence to some desired numerical accuracy. However, FPI generically suffers from stability issues: for a given point in parameter space, certain solutions may be unstable with respect to FPI, and so the algorithm will never converge on them. Here we review, following Ref. [18], the stability of FPI, first in the simple toy one-dimensional case, then generalising to the multi-dimensional case which is the case relevant for the CMSSM.
In one-dimension, FPI can be used to solve i.e. x * represents a fixed point of the function f (x). For FPI, we start with a guess x 1 of the parameter, and generate (hopefully better) successive approximations by Taylor expanding the right hand side around x * , we have, to first order in x * − x n , Stable convergence of FPI requires that, once x n is sufficiently close to x * , successive approximations are always closer: |x * − x n+1 | < |x * − x n |. Therefore, Conversely, if df (x * ) dx > 1, then FPI will be unstable around the solution. In our multi-dimensional case, x represents MSSM couplings and masses, and f (x) is the multi-dimensional function that represents what happens to x after running up to the high scale, applying boundary conditions, running back down and applying the other boundary conditions, i.e. one iteration: the loop shown in Fig. 1. In the following, bold capitals represent matrices, whereas bold font lower case letters are vectors. We now have N fixed point equations Proceeding similarly as we did in the one-dimensional case, we solve this system with successive approximations x n , where Taylor expanding the right hand side around x * , we have, to first order in (x − x n ), where D is the Jacobian matrix of f , evaluated at x * . Defining the matrix norm ||D|| of D, (3.8) as a stability condition. Conversely, from Eq. (3.7), we have If Eq. (3.10) is satisfied, the solution is clearly unstable to FPI 4 . Depending upon CMSSM parameters then, there may be solutions to the RGEs that have not been found by publicly available spectrum calculators, all of which rely on FPI. This was demonstrated to be true in Ref. [6], where one of the boundary conditions (Eq. (2.6)) was inverted: µ(M SUSY ) was fixed, and the equation was used instead to predict M Z . It was demonstrated that some points in CMSSM parameter space had several different values of µ(M SUSY ) that fit the empirically measured central value for M Z i.e. M Z (exp). These constituted multiple solutions of the RGEs for the same CMSSM point. Running FPI on some of the additional solutions showed that they were unstable, with parameters getting successively further from the solution in question. Table 1. Dimensionless order one input variables for each CMSSM point.

The Shooting Method
The shooting method [19] provides an alternative method to FPI of solving a multi-boundary problem. Instead of running the RGEs up and down, one runs only in one direction, by first applying the boundary conditions relevant for the starting scale and guessing any parameters that are not fixed by these boundary conditions. It is convenient for us to start at M GUT , applying the boundary conditions found in Eqs. We now face the fact that to truly have a chance of numerically finding all of the solutions, we must hunt in eleven dimensions for each CMSSM point. That is, for each choice of m 0 , M 1/2 , A 0 , tan β and sgn(µ), specifying a point in the CMSSM parameter space, there are 11 high-scale parameters (to be detailed below) which must be determined such that the low-scale constraints are satisfied. The overall strategy is to start at the high scale M GUT , input the MSSM parameters y i , consistent with M GUT -scale CMSSM boundary conditions, and evolve the RGEs down to M SUSY where the electroweak symmetry breaking boundary conditions should be met, then to the weak scale, where various weakscale boundary conditions should be met.
The MSSM renormalisation group equations may be phrased in general as where, i ∈ {1, 2, . . . , 109} and t = ln(Q/GeV), Q being the modified dimensional reduction scheme [15] We wish to find values of the V i such that the boundary conditions at M SUSY and M Z are satisfied: where the definitions of the f i are shown in Table 2. Thus, we must adjust the 11 V i Table 2. Output variables encoding the target boundary conditions. The various quantities are precisely defined in Appendix A. Note that f 11 measures how far we are away from a desired input value of tan β(input), whereas f 2 measures how far we are away from one of the Higgs potential minimisation conditions, i.e. Eq. (2.7).
in order to satisfy Eq. (3.12), providing a valid solution to the boundary conditions and RGEs (if such a solution exists). We can therefore view the system as 11 simultaneous nonlinear equations in the 11 V i . We solve them by two different methods, depending upon the context: Broyden's method [20] and a globally convergent Newton-Raphson method, which are both used to find roots of a multi-dimensional function. Here, we shall sketch the Newton-Raphson method that we employ; for more detail on it, or for information about Broyden's method, see Ref. [19]. We start with some guess for the initial parameters V i , and calculate the Jacobian matrix J ij = ∂f i /∂V j approximately by finite differences: 5 for small enough ∆V j . Since we solve j J ij ∆V j = −f i for ∆V j , which solves Eq. (3.14) to leading order for f i (V j + ∆V j ) = 0. Adding ∆V j onto V j should therefore take us closer to a solution of Eq. (3.12).
In practice, we require the step to be a descent direction, i.e. to first order is not satisfied in practice because we are too far from the minimum of f i f i and the second order terms in Eq. (3.14) give a sizeable correction, the algorithm backtracks until we find an acceptable step, as described in more detail in Ref. [19]. These updates are iterated until Eq. (3.12) is satisfied to a given numerical accuracy (specified here to be 10 −5 ). Unlike fixed point iteration, the shooting method, implemented using either Broyden or Newton-Raphson for the root finding, does not rely on any restrictions on the size of any derivatives in the vicinity of a solution. It is therefore able in principle to find a solution that is unstable to fixed point iteration. For given initial values of the V i , one solution (which has some 'catchment volume' in {V i } where the Broyden or Newton-Raphson method will converge to it, see Fig. 12) can be found. The idea then is to run with several randomly chosen V i in order to have the possibility of finding the usual fixed-point iteration solutions, the solutions found in Ref. [6] by scanning µ(M SUSY ), and any other solutions as well. We have the problem that, for a given number of initial V i , we can never guarantee that there are no other solutions other than those we have found. This is a general problem though in any numerical scan, and, failing some analytical breakthrough, there does not appear to be much we can do about it.

Numerical Procedure
We use a modified version of SOFTSUSY3.3.6 that implements the shooting method as described above. For a given CMSSM point, first we find a solution using the usual FPI method 6 . We then randomly perturb around this point and run the Broyden or Newton-Raphson root-finding algorithm in order to obtain a solution. The GUT scale parameters V i of this solution are checked against those corresponding to any solutions previously found. If any of the V i differs by at least 0.001, it is counted as an additional solution and saved. If the V i are all within 0.001, we save the solution with the smallest value of max i |f i |, since this represents the best numerical approximation to solving the boundary conditions. We repeat this procedure 2000 times with different random starting GUT scale initial parameters. It is hoped that 2000 times is sufficient to find all of the solutions, but we can never be sure that this is the case. Our results are usually phrased as scans across the CMSSM parameters, so statistical noise will be seen in any derived plots. Our random perturbations are based on draws from a Gaussian distribution centred on zero of width 1, which we denote G i (1) for each independent draw. Each of the V i are changed with a probability of 0.5. If perturbed, the V i are multiplied in turn by |1 + G i (1) × a i | (where we must pick values for the real coefficients a i ), with the exception of V 6 = [m 2 3 (M GUT )/10 6 GeV 2 ], which is instead multiplied by G i (1) × a 6 . V 6 is singled out for different treatment because we wish to consider either sign for its changed value.
The larger the a i , the less likely it is that SOFTSUSY will be able to return a physical solution (unphysical solutions could be due to the predicted presence of, for example, tachyons). On the other hand, larger a i explores a larger volume of MSSM parameter space within which extra solutions could be found. Choosing all a i of order one gives an unacceptably small efficiency, so as to make finding additional solutions extremely unlikely with 2000 'shots'. Indeed, trying to find any solution starting from a completely random starting point proved to be impossible: the efficiencies were much too small to find it in 2000 shots.
At tree-level, the M Z boundary conditions on the Yukawa and gauge couplings are specified unambiguously by Eqs. (2.1)-(2.4) and experimental data. The quantities µ(M SUSY ) and m 2 3 (M SUSY ) as derived from electroweak symmetry breaking can differ between multiple solutions, which in turn leads to electroweak one-loop corrections to sparticle threshold corrections to the gauge and Yukawa corrections. Some parameters a i correspond to boundary conditions that are fixed at tree-level, but may differ only in electroweak loop corrections. These a i we set to have a small value, whereas those a i corresponding to boundary conditions that may differ already at tree-level we set to have a larger value. By examining the efficiencies of finding solutions, we find, by trial and error, that the values a i = {0.03, 0.1, 0.03, 0.03, 1.0, 1.0, 0.03, 0.03, 0.03, 0.03, 0.1} give a reasonably high probability of finding the multiple solutions that we already know about, while retaining an efficiency (defined as the fraction of random perturbations resulting in a physical solution) that is not too small. The SOFTSUSY TOLERANCE parameter fixes how accurately the RGEs are solved. max i |f i | was also required to be less than TOLERANCE within 40 iterations of the Newton-Raphson algorithm (or Broyden's algorithm). We found that with higher values of TOLERANCE, we would obtain higher (and spurious) values of the number of multiple solutions at a CMSSM point: decreasing TOLERANCE further meant that some of the apparently different solutions were not good solutions after all, and some of them were actually inaccurately evaluated versions of the same solution. TOLERANCE was set to be 10 −5 for our results below. We have checked for many particular points with multiple solutions that reducing it further did not eliminate any such solutions (reducing TOLERANCE further resulted in a higher CPU time cost). We found that the statistical efficiency of our results did depend somewhat on whether we had implemented Broyden's algorithm, or Newton-Raphson, and furthermore on the precise implementation. Where there was a difference, we used the algorithm which found more multiple solutions, or which displayed less statistical noise in the results.

Regions of Multiple Solutions
We first provide a map of CMSSM multiple solutions, as determined by the multi-dimensional shooting method outlined above using Broyden's method. First, we show the results of a scan of m 0 and M 1/2 for A 0 = 0 and tan β = 40 (the plane investigated most thoroughly in Ref. [6]). We see that for for µ < 0, there are two disjoint regions with two solutions in Fig. 2a. In Fig. 2b, we see that for tan β = 40, A 0 = 0 and µ > 0 that there is a region at high m 0 with multiple solutions near the boundary of successful electroweak symmetry breaking. A line of three multiple solutions exists within this band. Remembering that our method is stochastic in nature, we are not surprised to see that this line is somewhat broken: we expect that, with more statistics, the gaps would be filled in (expecting that the multiple solutions are continuously connected in parameter space). Also, we see that to the top left hand edge of the multiple solution region, there are scattered points. We suspect too that here it is more difficult to find the multiple solutions, and that were we to have more statistics (i.e. run for a higher number of 'shots'), we would fill in this area. We also see multiple solutions near the electroweak symmetry breaking boundary at high m 0 for µ < 0 and tan β = 10 in Fig. 2c and for µ > 0 and tan β = 10 in Fig. 2d. In Fig. 2d, we display the 95% confidence level exclusion contour for ATLAS [21] and CMS [22] searches for events with hard jets and large missing transverse momentum in 5 fb −1 of √ s = 7 TeV LHC data. We see that this exclusion region covers much of the multiple solutions region.   Thus, as far as we can tell, exclusion limits derived upon those planes are not compromised by additional solutions.

Comparing algorithms
We now compare our results using the shooting method, to previous estimates in Ref. [6]. There, M Z was predicted by initially relaxing the boundary condition Eq. (2.6), scanning in one dimension, µ(M SUSY ), and using fixed point iteration in the remaining 10 dimensional parameter space per CMSSM point. Luckily, the existence of additional solutions was demonstrated, but in fact many of them could have been missed due to the fact that FPI was still being used, which could have the concomitant stability issues sketched in Sec. 3.1. In fact, the solutions found with the shooting method are a superset of those found in Ref. [6].
Comparing the results in Fig. 2 to equivalent estimates in Ref. [6] using the different methodology, we see some differences: we now find a much larger area of the parameter plane with multiple solutions. For example, we compare maps of a particular strip with two solutions in Fig. 3. We shall show below that this strip is of particular interest because it predicts sparticle masses that are not already excluded by LEP2. The new determination of this multiple solution region using the multi-dimensional shooting method is shown in green on the figure, whereas the one-dimensional method of Ref. [6] only finds the region in the black line. The insert shows that the width of the strip as measured by M 1/2 (for a given m 0 ) increases by a factor of several tens, providing a quantification of how much better the shooting method performs.
It is hard to compare the speed of the FPI computation versus that of our implementation of the shooting method, because often the shooting method will not find a solution. FPI also does not find additional solutions, and returns 'no-convergence' errors near the boundary of electroweak symmetry breaking. On the other hand when the shooting method works, the solution is typically much more accurate (i.e. M pred Z /M exp Z is typically much closer to 1). However, given these caveats, it is still probably helpful to give some idea of the relative speeds. In the case where FPI finds the same solution as our implementation of the shooting method, FPI is quicker. The vast majority of the run-time is taken up with running the RGEs between M GUT and M Z . Thus, we may use this is an approximate unit of time taken for the program to run. For each iteration of the shooting method we need to run the RGEs from M GUT to M Z 12 times (one for the central value, and 11 to determine the 11 forward derivatives in the Jacobian). This needs to be repeated typically a few times to converge on a solution (if indeed it does converge). FPI on the other hand, takes around 3-10 iterations, depending upon the parameter point 7 . Each of these iterations consists of running from M GUT to M Z and then back again, i.e. two time units. These rough estimates tell us that the shooting method is expected to be a few units slower than the FPI method when it achieves a solution. This is often increased by a factor of 2-10 because the shooting method takes that number of shots before a viable solution is found.
To get a quantitative understanding of the difference in speed of the algorithms, we performed rectangular scans over part of Fig. 6a for three algorithms: 1) pure FPI, 2) FPI to obtain a starting point, perturbing randomly around this point, then shooting with Newton-Raphson and a forward derivative, and 3) same as 2) but with a symmetric derivative (twice as many evaluations of the RGEs). We have tan β = 10, A 0 = 0 and look for positive µ solutions. Scanning with a 31 × 31 grid over the range m 0 = 1500 − 3000 GeV and M 1/2 = 500 − 800 GeV algorithms 1, 2 and 3 took 126, 437 and 701 seconds respectively. Compared with pure FPI, the shooting method will always take longer because it does all the computations of FPI (to find the initial guess), plus additional computation for Newton-Raphson convergence. It is thus best to summarise the running times as: the shooting method with forward derivatives added 247% to the running time of normal FPI, and with symmetric derivatives, 456%. The latter number is not quite twice the former (as naïvely expected) since, for example, the Newton-Raphson matrix inversion does not have to be repeated twice when going from forward to symmetric derivatives.
Because convergence properties differ over the CMSSM parameter space, we also scanned a rectangular grid that included part of the no-EWSB region. We ran a 31 × 21 grid over the range m 0 = 3000 − 4500 GeV and M 1/2 = 400 − 600 GeV (with all other parameters the same as in the previous scan). Per point, this region ran slower than the previous region.  pure FPI, these numbers correspond to an additional 162% and 311% for the forward and symmetric derivative shooting methods, respectively.

The LEP2 limit
We shall see that a sparticle mass bound derived at LEP2 in the context of the CMSSM imposes strong constraints upon the additional solutions. The 95% CL lower limit on the chargino mass from LEP2 is 102.5 GeV [24] within CMSSM-like models. The Lagrangian contains the chargino mass matrix −ψ −T Mψ +ψ + + h.c., whereψ + = for the charged higgsinos. We also have, at tree level, The chargino masses correspond to the singular values of Mψ + , i.e. the positive square roots of the eigenvalues of M †ψ + Mψ + [25]: This places a strong constraint on many of our multiple solutions, whose values of µ differ (and therefore the lightest chargino mass computed by Eq. (5.2) differs). In Fig. 4 we points that were already ruled out by the 95% CL region of ATLAS and CMS, as shown in Fig. 6, for example. The LHC limits are therefore unchanged on the tan β = 10, µ > 0 plane.
More generally, from Figs. 5, 6, 7 and 8 we see that every additional solution in the class adjacent to the no EWSB region fails the LEP2 chargino mass constraint for A 0 = 0, tan β = 10 and tan β = 40 and either sign of µ. This means that the only viable additional solution strip shown in these plots is the one at higher values of M 1/2 in Fig. 7b. We refer to this as the "phenomenologically plausible strip". We also show iso-contours of the lightest CP even Higgs mass in Figs. 5, 6, 7 and 8. One should bear in mind an estimated 3 GeV error upon the prediction coming from higher order corrections. Much of the available parameter space is ruled out by the recent LHC measurements of a boson consistent with a 125 GeV Higgs [26,27]. While this bound is not the main focus of the present paper, we see from Fig. 7b that there is an allowed region with multiple solutions: the phenomenologically plausible strip with m h 0 > 122 GeV.

The phenomenologically plausible strip
The only region of parameter space investigated that has additional solutions which are not excluded outright by the LEP2 chargino bounds is the thin phenomenologically plausible strip (the upper region in Fig. 7b, also shown in Fig. 3). In this region, the DR CP-odd neutral Higgs mass m A 0 tends toward zero, enhancing its loop correction in threshold effects (for example to Yukawa and gauge couplings). This enhancement allows for larger than usual corrections, leading to the presence of multiple solutions. For more detail see Ref. [6]. We shall now study the phenomenological properties of solutions along this strip in more depth. In particular, we wish to know how much the sparticle masses and other observables differ between the different solutions. If sparticle masses differ by less than 1% or so, then the solutions are similar enough such that LHC exclusion regions are unlikely to be different for the separate solutions. However, if there are significant differences in masses leading to different kinematics of SUSY signal events, cut efficiencies may be different enough to alter exclusion regions. We should bear in mind, however, that a future linear collider could be sensitive to fractional mass differences even at the sub-percent level, particularly in the slepton/chargino/neutralino sectors [28]. Thus, even sub-percent differences could be and the pseudo-scalar Higgs boson mass. We shall display some of these differences below. We shall also display branching ratios of stop and sbottom decays into charginos because they illustrate some significant differences between the two solutions. Fig. 9 shows how the fractional difference between the multiple solutions for various phenomenological properties varies across the strip for two fixed values of m 0 . The MSSM spectrum and coupling information was passed with the SUSY Les Houches Accord format [29] to SUSY-HIT [30], which calculated the branching ratios. Dark matter properties were then computed using micrOMEGAs [31,32]. For a given calculated phenomenological parameter X, the fractional difference is defined as |X 1 − X 2 |/X 1 where X 1 is the predicted value of X corresponding to the solution with a higher chargino mass and X 2 corresponds to the solution with the lower chargino mass. While the there are only sub-percent level differences between the masses of the lightest third-generation squarks and the lightest neutralinos, there are significant deviations for other quantities. The most significant differences are in the predicted thermal relic density of neutralino cold dark matter, Ω CDM h 2 , in which there can be 100% differences, depending upon the universal model parameters. There are also large differences at the tens-of-percent level in the branching ratio of the third generation squarks into charginos. At high m 0 , the chargino mass itself can also vary at the 10% level, which would necessarily affect the limits placed on this region of parameter space for analyses employing charginos in decays [33]. At lower values of m 0 , M χ ± 1 differs only at the percent level and so we would not expect there to be a large impact on limits. Thus, even though the chargino mass is large enough along the strip for the solutions to be viable, for a portion of this strip the phenomenological properties are so close that excluding one solution would also result in excluding the additional solution.
breaking at high m0, there is sensitivity in the lightest neutralino mass, but for lower m0, only the heavier neutralinos show any mass changes.

Example of Dark Matter Differences
Now we exemplify the most important difference between the solutions: that of a different predicted thermal relic density of dark matter in the universe. Recently, data from the Planck satellite have been used to derive the constraint [34] Ω CDM h 2 = 0.1198 ± 0.0026.
We place a dominant theoretical uncertainty on our prediction of 0.01 coming from loops (the thermal relic density is only calculated by micrOMEGAs to tree-level order), and therefore require the predicted thermal relic density of neutralinos to be Ω CDM h 2 ∈ [0.0998, 0.1398].
After a brief scan, we found a parameter point in the phenomenologically plausible strip (m 0 = 760 GeV, M 1/2 = 141.72 GeV, A 0 = 0, tan β = 40, µ < 0) where the standard solution predicts Ω CDM h 2 = 0.34, i.e. well outside of this range, but where the additional solution prediction of Ω CDM h 2 = 0.118 is near the central value. It turns out that this point has the χ 0 1 mass being approximately half of the Higgs mass. The χ 0 1 mass, which changes slightly between the solutions, is more exactly half the lightest CP even Higgs mass for the additional solution, which leads to very efficient annihilation of neutralinos through an s−channel h 0 into quark or lepton pairs, significantly reducing the relic density from 0.34 to 0.118.

Explorations in parameter space
It is possible to relax one or more of the constraints f i at the low scale and solve only a subset of these equations. Doing so will highlight the existence of multiple solutions, and give insight into how much the low-scale predictions change with respect to input parameters. For each constraint that is relaxed, there will be one high-scale parameter V j (it doesn't matter which one) that is no longer determined by the Newton-Raphson solver. Such high-scale parameter(s) are controlling parameters, in that we now fix them for the solution of that point. To get a feeling for how the parameter space looks we can scan these parameters -fix them at successive values -and plot the resulting function f i which has been left unconstrained.
Here, we choose to relax the M Z (pred) 2 constraint (i.e. f 1 can take any value) and the single controlling high-scale parameter is V 5 = [µ(M GUT )/1000 GeV] 2 . We scan V 5 over a large range and use the Newton-Raphson method to solve for the remaining 10 parameters using the remaining 10 constraints. Because the parameter space is large, we utilise a "line walking" technique in order to improve the efficiency of finding solutions for a given V 5 . Given two solutions that satisfy the 10 low-scale constraints, and are close in V 5 at the high-scale, we can use a linear approximation on the V i to project and make a guess at the high-scale parameters for the solutions either side of these original 2. Doing this successively allows us to "walk" through the high-scale parameter space, using V 5 as the parameter of the line, finding adjacent points that satisfy the 10 low-scale constraints. The initial point used is the one found by the standard FPI method.
The output of this line walking is a list of points in the 11-dimensional high-scale parameter space, one point for each value of V 5 . These points represent a function (f 1 =  M Z (pred) 2 /M Z (exp) 2 − 1) defined on the high-scale parameter space, and we plot all 11 projections of f 1 in Fig. 10, for both signs of µ (as illustrated by Fig. 10d 9 ). Points where the function f 1 crosses 0 correspond to full solutions of all boundary conditions and RGEs. Fig. 10d shows that, for this CMSSM point, there are two solutions for µ < 0 and one for µ > 0. These three solutions are evident in each of the other ten projections, although sometimes they are very close and hard to spot (for example in Fig. 10a). Each individual projection is non-smooth, and none of the one dimensional projections yield a uniquely defined function (i.e. for some values of each abscissa coordinate, there is more than Figure 11. Exploration of parameter space found by relaxing f 1 and f 11 in the CMSSM with for tan β = 40, A 0 = 0 GeV, m 0 = 2855 GeV and M 1/2 = 660 GeV. On the left we show a larger domain, but on the right we show a zoomed in domain of the µ < 0 region including two solutions of the full system. On the centre red contour line, M 2 Z is correctly predicted, on the centre blue contour, tan β(M Z (exp)) is. Where these contours intersect, all 11 constraints are satisfied; this occurs 3 times. The additional red contours on either side of the centre contour show where the error i.e. the predicted Z mass is 2 times the experimental one, physical or tachyonic, respectively. For the additional blue contours f 11 is ±5 × 10 −3 , i.e. tan β(M Z (exp)) is 0.5% higher or lower than the input value, as per Eq. (2.1). The white region has tachyonic A 0 s and does not yield physical solutions.
one predicted f 1 value). The non-smoothness yields practical difficulties when solving the system: the Newton-Raphson solver uses derivatives and the finite difference approximation to them becomes huge in the vicinity of a non-smooth piece. This can send the solver off to extreme and unphysical parts of parameter space, where we cannot obtain sensible solutions of the RGEs (for example, because of poles in renormalising quantities, or because tachyons are predicted). In this case, the attempt in question must be abandoned, and another started with a different random starting point. If that starting point is on the correct side of the kink, such that it can reach a full solution without encountering the kink again, then a solution may be found. This illustrates one of the reasons for requiring a stochastic approach where we fire many different shots in order to find the different solutions. To summarise: non-smoothness sometimes provide difficulties for the Newton-Raphson solver, whereas (as explained in Section 3.1), the magnitude of the derivatives at the horizontal zero line affects the stability of the FPI algorithm.
In a separate study we then relaxed two constraints: M Z (pred) 2 and tan β(M Z (exp)), equivalently f 1 and f 11 . We computed f 1 and f 11 over the two-dimensional plane of our freed-up parameters, chosen here to be M SUSY (M GUT ) and µ(M GUT ), controlled by V 11 and V 5 respectively. Within a rectangular region in this plane we pick random initial points and attempt to obtain a solution to the 9 low-scale constraints (we vary the density of points near the full solutions to provide more clarity). These functions each have contours at zero which predict the correct value for a low-scale constraint, in direct analogy with the lines in Fig. 10 passing through zero. Where these contours from the two unconstrained f 's intersect, we have a complete solution satisfying all 11 constraints. This is shown in Fig. 11. The standard solution for µ < 0 is at higher values of M SUSY (M GUT ). We can see from the figure that the additional solution here is close to the physical boundary, where the A 0 becomes tachyonic (shown as the white region in the plots). We also see from the figure that tan β(M Z (exp)) varies much more slowly across parameter space than M Z (pred) 2 . The figure also illustrates how difficult it is in general to obtain a solution: there is a large volume of unphysical (white) parameter space, and indeed the outer side of the red contours the Z 0 is tachyonic, yielding unphysical parameter space. The problem is much exacerbated in the full 11-dimensional case, rather than just this two-dimensional scan.
We also studied the behaviour of the Newton-Raphson solver by determining the 'catchment volumes' of the various solutions for one particular point in CMSSM parameter space. Figure 12 shows the result. We first find a solution using the FPI algorithm, and use this to fix 9 of the high-scale parameters. The other 2, being µ(M GUT ) and tan β(M GUT ), are scanned across a grid, as per the axis of the figure. Each point then corresponds to a distinct point in the high-scale parameter space, and we then use it as the initial guess for the Newton-Raphson solver. We then run the solver until it converges (or not) on a solution. If it converges, we colour the point depending on which solution it finds. In this example there are 4 different solutions, 2 with µ < 0 and 2 with µ > 0. The results are very sensitive to the particulars of the solver (such as tolerance, and criteria for convergence), but they serve to illustrate the fact that there are many sharp and non-trivial features in the equations that any solver must deal with.

Summary and conclusions
We have provided a method to search in an efficient, multi-dimensional way for multiple solutions to RGEs. The methodology that we employ should work in general for theories which have multiple renormalising quantities, with boundary conditions placed upon them at different renormalisation scales. It dispenses with fixed point iteration (the algorithm used by publicly available spectrum calculators), which can have stability problems for some of the solutions, and uses the shooting method instead. One shoots at low-scale boundary conditions only in one direction, by first imposing high scale boundary conditions and then running the RGEs down. Standard methods to solve non-linear simultaneous equations are then used -in this paper both Broyden's method and Newton-Raphson -to find solutions to the low-scale boundary conditions. To find several solutions, several shots are taken, each with random starting values of the high-scale parameters consistent with the high-scale boundary conditions. One problem that will be present in any such numerical scan is that one can never be completely sure that one has found all of the solutions, unless something analytically can be said about their number.
In a previous work [6] it was shown that the CMSSM has multiple solutions in some parts of its parameter space. The multiple solutions are characterised by the fact that they have the same m 0 , M 1/2 , A 0 , tan β and sgn(µ), yet are physically distinct, with different sparticle masses and couplings. The existence of multiple solutions was demonstrated by inverting one of the boundary conditions and using fixed point iteration in the other 10 dimensions of parameter space. Some solutions may be unstable to fixed point iteration, however, which means that such an algorithm will never find them. The shooting method has no such stability issues (that we are aware of), and can in principle find all of the solutions. This is done by randomly perturbing an initial guess of high-scale parameters consistent with the high-scale boundary conditions, then homing in on a solution with a non-linear multi-dimensional simultaneous equation solver. One repeats this several times, hoping with each 'shot' to find a different solution. Given enough shots, the method can find all of the solutions. However, for a finite number of shots, one can never guarantee that one has found all of them. In the present paper, we have shown that using the shooting method instead of fixed point iteration not only finds the previously discovered additional solutions, but also finds new ones too. It was also our intention here to study the multiple solutions phenomenologically, to see whether they are ruled out by experimental data, and to see to which extent they differ from the standard solutions (and if, for example, LHC searches need to be reinterpreted taking the multiple solutions into account).
We find that CMSSM multiple solutions come in two classes. One class, at extreme values of large m 0 just under the no electroweak symmetry breaking limit, is ruled out by LEP2 constraints on chargino masses. The other class is potentially phenomenologically viable, however. Studying this latter class, we see that sparticle masses are rather similar: most sparticles' masses are within 1% of the standard solution's. The most constraining LHC searches (various jets plus missing transverse momentum searches) are unlikely to have a large difference between solutions: if a CMSSM point is ruled out by such an analysis in the standard solution, it will be also ruled out for the additional solution (and vice versa). However, chargino masses may be around 10% different between solutions and so any CMSSM LHC analyses relying on charginos in decay chains are vulnerable to correction coming from the existence of the multiple solutions. Indeed, branching ratios of decays of stops and sbottoms into charginos are shown to differ by several tens of percent. Precision measurements at a future linear collider facility could certainly be sensitive to percent or per-mille differences in various sparticle masses, and so exclusion or measurement would need to take the additional solutions into account.
We predicted the thermal relic density of neutralinos, comparing the standard and the new, phenomenologically viable class of solutions. Here, we see large differences: a factor of 2-3 in Ω CDM h 2 is possible between the predictions of the two separate solutions. We have illustrated with a point that predicts far too much dark matter for the standard solution, where nevertheless the additional solution is near the observational central value derived by Planck. Thus, we conclude that analyses involving Ω CDM h 2 as a constraint should, of necessity, take multiple solutions into account. Many recent CMSSM analyses fit Ω CDM and other data (see, for example [35][36][37][38][39][40]), and these should be modified to include the additional solutions. It could be, for instance, that a better best-fit point exists within the additional solutions, meaning that the χ 2 fit frequentist type analyses' 95% confidence level contours move significantly.