Global fits of SUSY at future Higgs factories

In this work, we study the impact of electroweak and Higgs precision measurements at future electron-positron colliders on several typical supersymmetric models, including the Constrained Minimal Supersymmetric Standard Model (CMSSM), Non-Universal Higgs Mass generalisations (NUHM1, NUHM2), and the 7-dimensional Minimal Supersymmetric Standard Model (MSSM7). Using publicly-available data from the \textsf{GAMBIT} community, we post-process previous SUSY global fits with additional likelihoods to explore the discovery potential of Higgs factories, such as the Circular Electron Positron Collider (CEPC), the Future Circular Collider (FCC) and the International Linear Collider (ILC). We show that the currently allowed parameter space of these models will be further tested by future precision measurements. In particular, dark matter annihilation mechanisms may be distinguished by precise measurements of Higgs observables.


Introduction
After the discovery of the Higgs boson at the Large Hadron Collider (LHC), precisely measuring the properties of the Higgs boson is the next essential task for the high-energy physics community. For this purpose, various next-generation electron-positron colliders have been proposed, including the Circular Electron-Position Collider (CEPC) [1], the Compact Linear Collider (CLIC) [2], the Future Circular Collider (FCC-ee) [3], and the International Linear Collider (ILC) [4]. These state-of-the-art machines can not only scrutinize the nature of the Higgs boson but also shed complementary light on new physics, such as low-energy supersymmetry (SUSY; for recent overviews on the status of low-energy SUSY, see e.g. [5,6]). Restricted by collision energy, these colliders may not be able to directly produce supersymmetric particles, 1 but they can provide constraints on SUSY models through high-precision measurements of the Higgs and electroweak (EW) sectors.
As is well known, global fits can provide comprehensive information on new physics models, allowing us to infer the maximum amount of information on a given model from the widest range of experimental data [9,10]. Global fits assess and compare the validity of models, identify the ranges of model parameters with the highest likelihood or posterior probability and study the predictions and consequences for future searches. Consequently, global fits are an important part of the Technical Design Report (TDR) for future electron-positron colliders.
The GAMBIT Collaboration [11] have performed comprehensive global fits on five SUSY models: the Constrained Minimal Supersymmetric Standard Model (CMSSM) and its Non-Universal Higgs Mass generalisations (NUHM1 and NUHM2) [12], the sevendimensional weak scale phenomenological MSSM (MSSM7) with all parameters defined at the weak scale [13], and a four-dimensional electroweakino sector of the MSSM [14]. The likelihood functions of the first four global fits include several direct and indirect dark matter (DM) searches, a large collection of electroweak precision and flavour observables, direct searches for SUSY at the Large Electron-Positron collider (LEP), and Runs I and II of the LHC, and constraints from Higgs observables. All GAMBIT input files and generated likelihood samples for these models are publicly available online through Zenodo [15,16].
As a comprehensive global fit is computationally expensive, we take full advantage of the massive sets of publicly-available samples that GAMBIT generated by post-processing them with likelihoods of expected precision limits from future electron-positron colliders. By comparing the preferred regions and best fits before and after applying such likelihoods, we can estimate the prospective reaches of future colliders. This procedure is very different from previous global fits of SUSY with future electron-positron colliders (see, e.g., [17][18][19][20][21][22][23]).
The paper is structured as follows. We briefly introduce the theoretical framework of the CMSSM, NUHM1, NUHM2 and MSSM7 along with their parameters in Section 2. We then review the precision expectations of the CEPC, ILC and FCC-ee and present our post-processing strategy in Section 3. Each subsection in Section 4 contains implications for future collider searches on the global fit of one model. Finally, we summarize the results and draw our conclusions in Section 5.

The supersymmetric models
We restrict our consideration to specific scenarios within the CP-conserving, R-symmetric, minimal supersymmetric standard model (MSSM) [24]. The MSSM Lagrangian given in Sec. 5.4.3 of [11] defines the soft-breaking and SUSY-preserving parameters that describe the MSSM and fix our notation. Here we consider four distinct versions of the model where different constraints are applied to reduce the number of parameters. 2

CMSSM
The CMSSM is the most widely studied subspace of the general MSSM [25]. It is inspired by scenarios where SUSY breaking is transmitted through supergravity interactions, fixing the soft mass parameters at very high energy scales close to the Planck scale. Specifically, influenced by a minimal form of supergravity where the universal couplings are assumed, the universal scalar mass m 0 , the gaugino mass m 1/2 and the soft-scalar trilinear A 0 masses are defined at the Grand Unified Theory (GUT) scale, M X , through the following constraints: The Yukawa and gauge couplings of the MSSM can be fixed from data, while the Higgs sector parameters µ and b (often written as m 2 3 or Bµ) are partially constrained by the EW VEV, leaving the sign of the superpotential µ parameter and the ratio of the vacuum expectation values of the two Higgs doublets, as free parameters. In our work we input tan β at the scale m Z , following common conventions used in spectrum generators. This gives us four free continuous parameters and one free sign to fully specify scenarios in the CMSSM:

NUHM1
In this variant, the supergravity inspiration of the CMSSM is maintained but the strong assumption of minimality is relaxed a little. Specifically the GUT-scale constraint on the soft scalar Higgs masses is relaxed, introducing a new free parameter m H [26][27][28][29][30]. The constraint in Eq. 2 is replaced by The constraints in Eqs. 1, 3 and 4 are still applied in the NUHM1. This means the predictions of the NUHM1 are determined from five free continuous parameters and one sign,

NUHM2
Extending this idea slightly further, in the NUHM2 the constraint on the soft Higgs masses is even further relaxed so that m Hu (M X ) and m H d (M X ) become independent, real, dimension-one parameters at the GUT scale. Therefore, in the NUHM2 only the constraints in Eqs. 1, 3 and 4 are applied. This leaves six free continuous parameters and one sign to specify the physics predictions of the model,

MSSM7
In the MSSM7 the constraints are no longer applied at the GUT scale, allowing input parameters to be specified close to the weak scale and thus no longer requiring renormalization group (RG) running between this scale and the Planck scale to determine the models' physical predictions [31]. The number of parameters is reduced through the following constraints A e (Q) = 0 (14) where for the weak mixing angle we use a fixed value of sin , which is independent of precise extractions of the gauge couplings carried out during spectrum generation. Although in principle the scale Q can be arbitrary, Eqs. 11 and 12 are inspired by the impact of the GUT relation in Eq. 3 on the gaugino masses at the weak scale. So it does not make much practical sense to choose Q too far away from this scale, and in [32] and here we fix Q = 1 TeV. For the sfermion sector, we assume a common sfermion mass at 1 TeV with no flavour changing and assume only the third generation up-and down-type soft trilinear contribute non-negligibly. In the MSSM7 we also make the trade between VEVs and |µ|, b. As in the NUHM2 both soft Higgs masses are free inputs, m Hu (Q) and m H d (Q), albeit with the parameters input at very different scales which will substantially change the prior distributions they are drawn from. Similarly the sfermion masses are also no longer split by the RG flow, and the soft trilinears are treated differently with A t and A b as free parameters and A τ = 0. The final list of free parameters then has seven continuous parameters and one sign: Global fits of these four models have previously been performed with GAMBIT in [12,13]. Here we will post-process the results from these global fits to explore the impact of future electron-positron colliders. For the samples we use, the constraints listed above have been applied in FlexibleSUSY [32,33], which calculates the pole masses and DR couplings used in the global fit. The values of the parameters at different scales are connected by two-loop RGEs which FlexibleSUSY obtains 3 with the help of SARAH [36][37][38][39]. The scanned parameter ranges can be found in [12,13]. All dimensionful parameters are allowed to vary up to 10 TeV, safely covering scenarios that are motivated by the hierarchy problem and the most phenomenologically interesting regions within reach of colliders.

Study Strategy
In our analysis, we post-process the publicly available data on Zenodo for global fits of the GUT scale SUSY models [15], and MSSM7 [16], with additional likelihoods for precision measurements of the SM Higgs observables at the proposed Higgs factories. The total likelihood is thus, L = L GAMBIT · L Higgs factories (16) where the likelihood already computed in the publicly available data, L GAMBIT , includes contributions from a large collection of present constraints on dark matter, electroweak precision, flavour observables, sparticles and the SM-like Higgs boson. Some of the constraints have been improved since the publication of the data, but they will not qualitatively affect our calculation here. Since we post-process the data, the parameter ranges and priors are the same as in the original studies [12,13]. A Higgs factory with e + e − collisions at a center-of-mass energy of 240-250 GeV exploits the Higgsstrahlung process e + e − → hZ.  Table 1: Estimated statistical precision for Higgs measurements at the proposed ILC program with various center-of-mass energies [40], the FCC-ee program with 5 ab −1 integrated luminosity at √ s = 240 GeV and 1.5 ab −1 integrated luminosity at √ s = 365 GeV [41,42], and the CEPC program with 20 ab −1 integrated luminosity at √ s = 240 GeV and 1 ab −1 integrated luminosity at √ s = 360 GeV [43].
With millions of Higgs bosons produced and the clean experimental conditions at lepton colliders, both the inclusive cross section σ hZ independent of the Higgs decays, and the exclusive channels of individual Higgs decays in terms of σ hZ × Br, can be measured to impressive precision. The cross sections of vector boson fusion processes for the Higgs production are relatively small at low values of center-of-mass energy. Only the main decay modes can be measured. With the center of mass energy increasing, the cross sections of vector boson fusion processes grow logarithmically and can provide crucial complementary information.
In Table 1, we list the anticipated precision of measurements of the Higgs boson rates in the running scenarios of various machines in terms of their center of mass energies and the corresponding integrated luminosities. For the center-of-mass energy of 240−250 GeV, only the measurement of h → bb is provided in the vector boson fusion process, which is listed in the last row.
We define the new likelihoods for the proposed Higgs factories simply as where (19) and the index i runs over all the Higgs search channels in Table 1. We take the values of branching ratios for the 125 GeV SM Higgs as Br SM i , listed in Table 2. The total uncertainties σ m h and σ µ i consist of experimental uncertainties and theoretical uncertainties, σ tot = σ 2 the + σ 2 exp . The parametric uncertainties are involved in through nuisance input parameters. For the SM Higgs mass, the present experimental uncertainty, σ exp m h = 0.17 GeV [44], is already small and can be neglected compared with the theoretical uncertainty in GAMBIT, σ the m h = 2 GeV. Thus, we ignore the contribution Decay mode Branching ratio Theoretical error Table 2: SM predictions of the decay branching ratios for the 125 GeV Higgs boson [50], as well as the relative theoretical uncertainties.
from the Higgs mass measurement at future Higgs factories to avoid double counting the likelihood for the Higgs mass. We use the anticipated precisions displayed in Table 1 as the experimental uncertainties on future signal strength measurements, while the theoretical uncertainty is in the same order and therefore cannot be ignored. In Table 2, we show the present theoretical uncertainties of branching ratios for the 125 GeV SM Higgs boson, which are expected to be improved in the future. We consider theoretical uncertainties σ the µ i = kσ SM µ i for BSM predictions and by default we use k = 0.2 in our likelihood, representing a scenario where the theoretical uncertainties are reduced to the same level as the expected experimental uncertainties. However, the uncertainties on the branching ratios computed for BSM models in GAMBIT, which use SUSY-HIT 1.5 [45][46][47] via DecayBit [48], are obviously larger than those in the SM, as carefully shown in [49] and discussed in [23]. We discuss the impacts of SUSY contributions to the uncertainties, and using different choices of k in Section 4.2.
As for the central values of future signal strength measurements, µ obs i , in our scenario we assume that signal strengths predicted by the best-fit (BF) points in the GAMBIT datasets are measured. Since the GAMBIT data are sampled with Diver [51] and Multi-Nest [52,53], the regions around the BF points are densely sampled, such that our postprocessing procedure should be reasonable. Since the BF point agrees exactly with our assumed measurements, the BF point should not change through post-processing. The extent to which the confidence regions around the BF shrink, however, will reveal the potential impact of precise Higgs factory measurements. We will also use other central values to check the validity of conclusions drawn from this choice in Section 4.1.
The proposed electron-positron colliders are also designed to run at the Z-pole, and have an excellent capability for precise measurements of EW observables. These measurements are complementary to the Higgs boson coupling measurements. For instance, in the so-called "blind spot" the coupling between the lightest stop and the light Higgs boson vanishes, but the stop contribution to the EW precision observables can be visible [19]. In the existing GAMBIT data, the strong coupling at the scale m Z in the M S scheme α MS s (m Z ) and the top pole mass m t are nuisance parameters varying within ±3σ of their observed central values, while the Z pole mass is fixed to 91.1876 GeV. Beside, the W mass and sin 2 θ W are outputs calculated in PrecisionBit [48]. Therefore, we also   [54], which are originally from [55][56][57][58].
investigate the impact of the expected EW precision measurements using For the same reason as the L Higgs factories , we set the central values of these measurements to those predicted by the BF point. They are listed for the CMSSM in Table 3, along with the central values used in the likelihood of previous GAMBIT fits and the anticipated precisions of measurements of the EW observables at future lepton colliders. As for theoretical uncertainties, we adopt 1.5 × 10 −5 for sin 2 θ W and 1 MeV for m W .
In total, we post-processed 7.1 × 10 7 viable samples for the CMSSM, 9.4 × 10 7 samples for the NUHM1, 1.2×10 8 samples for the NUHM2, and 1.8×10 8 samples for the MSSM7. For each of the models the post-processing took several days on 1280 supercomputer cores. In order to compare our results with those shown in previous GAMBIT papers [12,13], we present and plot our results in the same way. The samples are sorted into 60 bins across the range of data values in each dimension, and the resulting profile likelihoods are interpolated with a bilinear scheme in Pippi 2.0 [59].
The post-processing is performed using GAMBIT 2.1.1, while the original global fits employed GAMBIT 1.0. The spectrum generators in the two versions are slightly different, as GAMBIT 2.1.1 includes bug fixes in calculations of the Higgs branching ratios. Therefore, the branching ratios of the SM-like Higgs in the two results are not fully consistent. Fortunately, the difference is slight and can be neglected compared to the current precision of Higgs measurements at the LHC. Using the new branching ratios, the best-fit regions in the old studies will not change qualitatively. Thus, it is still reasonable to use the results shown in papers [12,13] as the global fits of the present likelihood.

Results
In this section, we compare profile likelihoods with and without the additional likelihood for the Higgs measurements at future electron-positron colliders, by taking CEPC as an example, in the CMSSM, NUHM1, NUHM2 and MSSM7. Moreover, the dependence of the results on assumptions about central values of Higgs measurements at future facilities and theoretical uncertainties are investigated. We also compare the sensitivities of the CEPC, FCC-ee and ILC.  The profile likelihood ratio in the CMSSM for the present constraints (left panels, taken from [12]) and including additional CEPC measurements (right panels), with 68% and 95% CL contours drawn in white, and the best-fit point indicated by a star. We use the branching ratios of the best-fit point in the CMSSM as the central values of the measurements at the CEPC, and the theoretical uncertainties are k = 0.2 times smaller than the current SM Higgs theoretical uncertainties.

CMSSM
We show the 2D profile likelihoods in Figure 1 for the input parameters of the CMSSM without (left panels) and with (right panel) the additional likelihood for the Higgs cou-pling measurements at the CEPC. Here we assume that the central values of measurements at the CEPC are same as those predicted by the best-fit point of the CMSSM, and the theoretical uncertainties are k = 0.2 times smaller than the current SM Higgs theoretical uncertainties.
We see that a large part of the region favored by present constraints is excluded by the Higgs precision measurements at the CEPC. As we use the best-fit point (the white star in the plots) to set the central values of the measurements at the CEPC, it leads to −2 ln L CEPC = 0 at the best-fit point, so the position of the best-fit point holds still, and the preferred regions shrink significantly towards the best-fit point. It is very useful to classify the regions according to the possible DM annihilation mechanisms as follows: • stop co-annihilation: where 'heavy' means a heavy Higgs like A 0 or H 0 . The best-fit point is located in the stop co-annihilation region. The stop co-annihilation region prefers large and negative A 0 , and extends below the lower bound of bottom panels of Figure 1. However, vacuum stability problem need to be examined carefully in this region. We see that the regions of large m 0 , m 1/2 and tan β, i.e. the A/H-funnel and chargino co-annihilation regions vanish. Besides, the sign of µ in the remaining stop co-annihilation regions is always negative.
The sign of the µ parameter highly affects several physical observables, such as the recently updated anomalous magnetic moment of the muon a µ [60]. As mentioned in [12], the fit favors µ < 0 versus µ > 0 by ∆ ln L = 0.4, mainly because of the LHC Higgs signal likelihood. With the impressive precision of Higgs property measurements at the CEPC, the distinction between µ > 0 and µ < 0 reaches more than 2σ.
Note that the classification of dark matter annihilation mechanisms is not exclusive. A sample can lie in more than one region. The classifications provide information about the relationship between sparticle masses. Samples in the same region have similar parameter values and mass spectra, as well as similar SM Higgs branching ratios.
In order to further understand the significant impact of the CEPC, as shown in Figure 1, we display the SM-like Higgs decay branching ratios predicted by the best-fit points in each of the regions in Table 4. The combined experimental and theoretical uncertainty is also shown. The parameter values, mass spectra, and present likelihood contributions of these best-fit points can be found in [12]. We see that the differences in BR(h → bb), BR(h → W W * ) and BR(h → ZZ * ) between the best-fit points of the stop co-annihilation region and the A/H-funnel region or theχ ± 1 co-annihilation regions are significantly larger than the corresponding total uncertainties. This is the reason why the A/H-funnel region and theχ ± 1 co-annihilation region are excluded when we assume that the CEPC measures exactly the central values predicted by the best-fit point in the stop co-annihilation region.
It is obvious that the results shown in Figure 1 and the above conclusions depend on assumptions about central values of the Higgs measurements at the CEPC. Therefore, in Figure 2 we display the 2D profile likelihoods assuming the central values of the  Table 4: SM-like Higgs signal strengths and normalized cross sections of the best-fit points for the present likelihood in the CMSSM, for each of the regions characterised by a specific mechanism for suppressing the relic density of dark matter. We also give the total uncertainty used in L CEPC in the last column, assuming k = 0.2.
Higgs measurements at the CEPC to be the values of the best-fit point in each DM annihilation region, not just the overall best-fit point. Whereas the A/H-funnel region predicts SM-like Higgs couplings, theχ ± 1 and stop co-annihilation predictions for the CEPC are expected to be in about 1σ and 5σ tension with the SM, respectively.
It can be seen that the favored regions change dramatically with the central values changing from the values of one best-fit point to another. As expected, the center of the favored regions is on the chosen best-fit point. In the bottom left panel of Figure 2, where the central values are the same as the best-fit point in the A/H-funnel region, the favored regions are not as narrow as before. There is a small stop co-annihilation region inside the 95% CL region. The samples in the 95% CL region of m 0 < 8.5 TeV and m 1/2 > 2.6 TeV mainly satisfy theχ ± 1 co-annihilation condition, while the 68% CL region is almost pure A/H-funnel region. In the whole region, the µ parameter is always negative. The confidence regions found when setting the central values of the CEPC measurements to those predicted by the best-fit point in theχ ± 1 co-annihilation region, shown in the bottom right panel of Figure 2, are even wider. All three DM annihilation regions exist in the 68% CL region. However, there are no samples with µ > 0 in the 95% CL region.
The sfermions in the A/H-funnel and theχ ± 1 co-annihilation regions are heavier than about 5 TeV. For instance, the best-fit point in theχ ± 1 co-annihilation region has mt 1 mτ 1 6.5 TeV and mq 10 TeV for the first two generation sfermions. The supersymmetric contributions to the SM-like Higgs branching ratios tend to decouple in the high mass region. Nevertheless, the SM-like Higgs branching ratios in these regions remain quite different from the SM values listed in Table 2. Various uncertainties and  Table 2. To some extent, the SM-like Higgs branching ratios of a point in the high mass region can be treated as the SM values, such as the best-fit point of the A/H-funnel region or theχ ± 1 co-annihilation region. We also show the 2D profile likelihood for EW measurements at the CEPC in Figure 3, where the left panel implements the EW measurements on their own and the right panel combines likelihoods for the EW and Higgs measurements. The EW measurements visibly narrow the favored regions, though not to the same degree as the Higgs measurements. The two kinds of measurements complement each other well, as combining them gives extremely strong constraints on the parameter space of the CMSSM, shown in the right panel of Figure 3. The BF point's predictions for EW observables are expected to be in about 3σ tension with the SM.

NUHM1 and NUHM2
In the right panels of Figure 4 we show the joint profile likelihood ratio including the proposed CEPC constraints for the input model parameters of the NUHM1 and NUHM2, accompanied by the profile likelihood ratio without L CEPC in the middle panels and the dark matter annihilation mechanisms in the left panels. The definitions of the dark matter annihilation regions are same as those in the above subsection, except for an additional category • stau co-annihilation: mτ 1 ≤ 1.2 mχ0 1 , which is absent in the CMSSM results. With the extra freedom presented in the Higgs sector of the Non-Universal Higgs Mass models, the µ parameters decouple from m 0 , leading to arbitrarily light Higgsinos. Thus, the chargino co-annihilation region expands significantly. Meanwhile, the best-fit points in both the NUHM1 and NUHM2 results are also located in the stop co-annihilation region and have slightly larger likelihoods than the best-fit point in the CMSSM.
As in the CMSSM, with the addition of L CEPC , the preferred regions shrink significantly towards the best-fit points in the NUHM1. However, the stop co-annihilation region here overlaps with all other three regions in all of the parameter planes. As a result, the remaining 2σ regions also contain some chargino co-annihilation region. Note that the 1σ region is a pure stop co-annihilation region. Samples of µ > 0 are all excluded. In comparison to the NUHM1 results, the NUHM2 results show larger 1σ regions but smaller 2σ regions, with no chargino co-annihilation region inside the 2σ region.  Figure 4: Profile likelihood ratio, without the CEPC likelihood (middle panels, taken from [12]) and with the CEPC likelihood (right panels), for the NUHM1 (top two rows) and the NUHM2 (bottom two rows). Colour-coding in the left panels (taken from [12]) shows the mechanisms active in models within the 95% CL contour for avoiding thermal overproduction of neutralino dark matter. The overall best-fit point is indicated by a white star, while the best-fit points in each region are indicated by colored stars. The assumptions about the CEPC likelihood are the same as those in Figure 1.
The extent of the confidence regions surrounding the best-fit point depends to some degree on the assumptions made about the theoretical uncertainties on the Higgs branching ratios. Comparing with the top right panel of Figure 4 which uses 0.2 of the current theoretical uncertainties in L CEPC , we show the result without theoretical uncertainties and with full current SM theoretical uncertainties in the top left and top right panel of Figure 5, respectively. In general the confidence region shrinks a lot with smaller theoretical uncertainties, though the dark matter annihilation mechanisms and the sign of µ in the 1σ regions remain the same. Without theoretical uncertainties, only small regions around the best-fit point remain. Assuming no improvement on the theoretical uncertainties, the favored regions still shrink, but this would substantially negate the advantage of Higgs factories. In addition, we add SUSY contributions to the theoretical uncertainties according to Ref. [23], and display the results in bottom panels of Figure 5.  In Figure 6, we show the results on the mχ± 1 -mχ0 1 plane. We see that as the dark matter annihilation mechanisms constrain the dark matter mass and the relationships between sparticle masses, and that the masses of relevant sparticles are further restricted into limited ranges by L CEPC . In the stop co-annihilation region, the wino-dominated chargino mass is about twice as large as the bino-dominated dark matter mass, because

MSSM7
The results of the MSSM7 scan are shown in Figure 7, on the parameter planes of µ-M 1 , M 2 -mf and µ-tan β, and on the mass plane ofχ 0 1 -t 1 . Here µ and M 1 are presented at the scale M SUSY = √ mt 1 mt 2 . Two new regions • sbottom co-annihilation: mb where 'light' may be h or Z. In the sbottom co-annihilation region, the lightest sfermion ist 1 and therefore it highly overlaps with the stop co-annihilation region. In contrast to the GUT-scale models, the overall best-fit point of the MSSM7 is located in the chargino co-annihilation region, and the corresponding best-fit likelihood is improved.
We see that the 1σ ranges are shrunk significantly, similar to the results of GUTscale models, but the 2σ ranges of parameters are not visibly reduced by the CEPC constraints, except for the µ parameter. With free m 2 Hu and m 2 H d , the lightest chargino and neutralino can be higgsino-dominated for almost any value of M 2 , mf , A d 3 , A u 3 and tan β. Considering that we set the central values of Higgs measurements at CEPC to the  Figure 7: Profile likelihood ratio without the CEPC likelihood (middle panels, taken from [13]) and with the CEPC likelihood (right panels) for the MSSM7. Colour-coding in the left panels (taken from [13]) shows the mechanisms active in models within the 95% CL contour for avoiding thermal overproduction of neutralino dark matter. The overall best-fit point is indicated by a white star, while the best-fit points in each region are indicated by colored stars. The assumptions about the CEPC likelihood are the same as those in Figure 1.
values of the overall best-fit point, most of the chargino co-annihilation region survives. Thus, the ranges of those parameters cannot be restricted. On the other hand, the µ parameter has to be lighter than M 1 in the chargino co-annihilation region, so that the lightest chargino is not bino-dominated, shown in the top left panel of Figure 7. When the other regions are excluded, the upper limit on µ drops accordingly. As the A/H funnel region and the chargino co-annihilation region overlap heavily in all the planes, part of the A/H funnel regions escapes from the restriction of the CEPC Higgs measurements. Besides, there are a few samples satisfying the stop and sbottom co-annihilation conditions, which can be seen in the low panels of Figure 7, but they fulfill the A/H funnel selection at the same time. We also see one sample satisfying the h/Z funnel and the chargino co-annihilation conditions simultaneously. Overall, Higgs measurements at the CEPC do not have much power to discriminate between different DM annihilation mechanisms in the MSSM7.
We present a brief comparison for the potential reach of the CEPC, FCC-ee and ILC in Figure 8, together with the top left of Figure 7. The results include measurements at all proposed center-of-mass energies listed in Table 1. The CEPC and ILC result in slightly stronger constraints than the FCC-ee and almost same contour regions, but for different reasons. The CEPC proposes a higher integrated luminosity at √ s = 250 GeV, while the ILC has more center of mass energy options. The difference of favored regions between these facilities are small due to the theoretical uncertainties which are larger than the anticipated precision in some signal channels.
Before concluding, we consider the implications on our models of two outstanding experimental anomalies. First, the surviving regions shown for all of these SUSY models cannot make a significant contribution to resolving the recently updated muon g − 2 anomaly [60], as the unified sfermion mass mf is pushed to a relatively large value by other constraints, such as LHC sparticle searches and B-physics constraints. An MSSM explanation of the muon g − 2 anomaly at the 2σ level requires light EWinos and light sleptons (see, e.g., [61][62][63][64]), which may be accessible at the future runs of the LHC [64,65].
Second, the new CDF II measurement of the W -boson mass shows a 7σ deviation from the SM prediction [66]. This differs significantly from the central value used for m W in previous GAMBIT fits. Therefore, all samples in our favored regions disagree with the new measurement by at least 2σ. In any case, it is difficult to explain the large deviation in the general MSSM [67][68][69].
Finally, let us remark on the issue of naturalness in our scenarios. We quantified the fine-tuning for the best fit points in the GUT-scale SUSY models through Using SuSpect [70], we found about 4000 for the CMSSM, 700 for the NUHM1 and 900 for the NUHM2. The fine-tuning in the favored regions of the CMSSM is generally larger than that in the NUHM1/2, because the condition m 2 0 = m 2 Hu = m 2 H d leads to a strict constraint on m 0 . Investigating the fine-tuning cost comprehensively is beyond the scope of this paper, but the fine-tuning of these benchmarks should be a reasonable indication of the typical level of tuning one expects in the surviving samples.

Conclusions
Utilizing the publicly available data for SUSY global fits from the GAMBIT community, we examined the potential impact of measurements at proposed Higgs factories on several constrained versions of the MSSM, namely the CMSSM, NUHM1, NUHM2 and MSSM7. We post-processed all the samples to calculate the likelihood for Higgs measurements from the CEPC as an example, and then compared profile likelihood ratios with and without this additional likelihood from the CEPC. The preferred regions in these models are significantly shrunk by the precise Higgs measurements at the CEPC. As a result, the possible dark matter annihilation mechanisms in the models and signs of µ parameter could be distinguished by CEPC measurements. Comparing results in different models, the constraints on model parameters are weaker in models with more input parameters, i.e., looser correlations between model parameters. The specific favored and excluded parameter regions highly depends on the assumed central values of the Higgs measurements, and are mildly dependent on assumptions about theoretical uncertainties. Since the projected experimental uncertainties are even better than the current theoretical uncertainties, reaches of Higgs factories are dominated by the theoretical uncertainties. Reducing these theory uncertainties is an important and significant challenge. Under the assumption that the future theoretical uncertainties could be reduced to the same level of expected experimental uncertainties, we compared the impact of the CEPC, FCC-ee and ILC on the MSSM7. We found that their reaches are quite similar, and slightly better with higher proposed luminosity or higher center-of-mass energy, as expected. In summary, with high-precision Higgs coupling measurements, future Higgs factories can significantly advance our understanding of the MSSM parameter space and mass spectrum, and could be complementary to dark matter searches and EW precision measurements.