Daniele Barducci
email: barducci@lapth.cnrs.fr

Alexander Belyaev
email: a.belyaev@soton.ac.uk

Aoife K M Bharucha
email: aoife.bharucha@cpt.univ-mrs.fr

Werner Porod
email: porod@physik.uni-wuerzburg.de

Veronica Sanz
email: v.sanz@sussex.ac.uk

Uncovering Natural Supersymmetry via the interplay between the LHC and Direct Dark Matter Detection

come   L'archive ouverte pluridisciplinaire

Contents 1 Introduction

The naturalness of Supersymmetry (SUSY), which has been subjected to much discussion and thorough investigation for more than two decades , has become even more relevant now, as the Large Hadron Collider (LHC) collaborations ATLAS and CMS have started to probe SUSY in the TeV region. Indeed, the lack of evidence for superparticles at the CERN LHC, along with the rather high value of the Higgs boson mass in the context of the Minimal Supersymmetric Standard Model (MSSM), raises the question of whether the remaining allowed parameter space suffers from a high degree of fine-tuning, and if there is any parameter space of Natural SUSY (NSUSY) left. We discuss this problem in the framework of the well motivated MSSM.

Based on standard measures of fine tuning [START_REF] Ellis | Observables in Low-Energy Superstring Models[END_REF][START_REF] Barbieri | Upper Bounds on Supersymmetric Particle Masses[END_REF], the NSUSY parameter space was originally associated with light higgsinos, gluinos and stops, the SUSY partners of the SM Higgs, gluons and top quark. The present LHC limits on the masses of the latter two are approaching the TeV scale, under the assumption that the mass gap between these states and the lightest supersymmetric particle (LSP) , χ0

1 , is large enough (see for example [START_REF] Chatrchyan | Search for gluino mediated bottom-and top-squark production in multijet final states in pp collisions at 8 TeV[END_REF][START_REF] Chatrchyan | Search for supersymmetry in hadronic final states with missing transverse energy using the variables AlphaT and b-quark multiplicity in pp collisions at 8 TeV[END_REF][START_REF] Aad | Search for supersymmetry at √ s=8 TeV in final states with jets and two same-sign leptons or three leptons with the ATLAS detector[END_REF][START_REF] Aad | Search for squarks and gluinos with the ATLAS detector in final states with jets and missing transverse momentum using √ s = 8 TeV proton-proton collision data[END_REF][START_REF] Aad | Search for top squark pair production in final states with one isolated lepton, jets, and missing transverse momentum in √ s =8 TeV pp collisions with the ATLAS detector[END_REF]). Here χ0

1 is the lightest of the four neutralinos, the mass eigenstates arising from the mixing of the fermionic component of the Higgs and gauge superfields, ( H0 d , H0 u ) and ( B0 , W 0 ), which are commonly known as the higgsinos, the bino and the wino. Moreover, these limits have a certain degree of model dependence and for stops, t1,2 , the experimental limits rely on certain decay channels being dominant (e.g. t1 → t χ0

1 , with t the top-quark) or on a substantial mass splitting between t1 and χ0

1 , and can be significantly relaxed. For example, in the scenario under consideration in this paper with higgsino like dark matter (DM), the stop branching ratios strongly depend on the left-right admixture of the lightest stop. Therefore model independent collider bounds on stops are weak or non-existent. It has however been shown that usual fine-tuning measures, defined as the sensitivity of the weak scale to fractional variations in the fundamental parameters of the theory, can be low even if the masses of the supersymmetric scalars are large. This happens in the so called "hyperbolic branch"(HB) [START_REF] Chan | Naturalness, weak scale supersymmetry and the prospect for the observation of supersymmetry at the Tevatron and at the CERN LHC[END_REF] or "focus point" (FP) [START_REF] Feng | Focus points and naturalness in supersymmetry[END_REF][START_REF] Feng | Multi -TeV scalars are natural in minimal supergravity[END_REF][START_REF] Feng | Focus Point Supersymmetry Redux[END_REF] regions of the minimal super gravity (mSUGRA) parameter space, where the value of the Higgs mass parameter, µ, can be low if the universal gaugino mass M 1/2 is not too large, as a consequence of the subtle interplay between the electroweak (EW) gauge couplings and the top-Yukawa coupling in the evolution of the squared Higgs mass parameters using the renormalization group equations (RGE). Moreover it was recently argued [START_REF] Baer | How conventional measures overestimate electroweak fine-tuning in supersymmetric theory[END_REF] that EW fine-tuning in SUSY scenarios can be grossly overestimated by neglecting additional terms, stemming from the ultra-violet (UV) completion of the model, that can lead to large cancellations favouring a low µ parameter, but not necessarily a low stop mass up to a certain limit. Taking this point of view, we will take a low µ parameter to be the definition of NSUSY throughout our study.

In the case µ M 1 , M 2 (the EW gaugino mass parameters) one finds that the three lightest neutralino and chargino mass eigenstates, χ0

1 , χ0 2 and χ± 1 , are quasi-degenerate and that these states are nearly pure higgsinos. In this scenario the DM relic density is typically below the WMAP [START_REF] Hinshaw | Nine-Year Wilkinson Microwave Anisotropy Probe[END_REF] and PLANCK [START_REF] Ade | Planck 2013 results. XVI. Cosmological parameters[END_REF] measurements, because of the high rate of higgsino annihilation to standard model (SM) gauge and Higgs bosons and the higgsino coannihilation processes [START_REF] Feng | Prospects for indirect detection of neutralino dark matter[END_REF][START_REF] Baer | Neutralino relic density in minimal supergravity with coannihilations[END_REF]. One should note that this parameter space with the relic abundance below the experimental constraints is not however excluded, since the remaining relic abundance can be accounted for by other additional sources, e.g. axions. This NSUSY scenario, which is characterised by relatively light higgsinos in comparison to other SUSY particles, is not just motivated by its simplicity, but also by the lack of evidence for SUSY to date. We take advantage of the fact that NSUSY scenarios can be effectively described by a two dimensional parameter space, defined by the DM mass, i.e. the mass of χ0

1 , and ∆M, the mass difference between the DM candidate and the next to lightest supersymmetric (NLSP), typically χ± 1 , and our study explores the complementarity of the LHC and direct detection (DD) DM search experiments in covering this region. Such complementarity was the subject of recent studies, see e.g. [START_REF] Crivellin | Light stops, blind spots, and isospin violation in the MSSM[END_REF] and references there in.

It was already shown a decade ago that the HB/FP parameter space is challenging to probe at the LHC [START_REF] Baer | Indirect, direct and collider detection of neutralino dark matter[END_REF] even if the mass gap between χ0

1,2 and χ± 1 is large enough to provide leptonic signatures. The most challenging case arises when the mass gap between these states is too small to produce any detectable leptons. The only way to probe such a scenario is via mono-object signatures, i.e. signatures involving a high transverse momentum particle recoiling against missing transverse emerge (E miss T ), of which the monojet signature is particularly of relevance at the LHC, as initially suggested in [START_REF] Alves | Where the Sidewalk Ends: Jets and Missing Energy Search Strategies for the 7 TeV LHC[END_REF] for generic compressed spectra.

This technique has already been used in studies of quasi-degenerate higgsino spectra via monojet+E miss T and monojet+E miss T +soft di-lepton signatures for the NSUSY parameter space we consider [START_REF] Han | Probing Light Higgsinos in Natural SUSY from Monojet Signals at the LHC[END_REF][START_REF] Schwaller | Compressed electroweakino spectra at the LHC[END_REF][START_REF] Han | Hunting Quasi-Degenerate Higgsinos[END_REF][START_REF] Baer | Monojets and mono-photons from light higgsino pair production at LHC14[END_REF][START_REF] Baer | Monojet plus soft dilepton signal from light higgsino pair production at LHC14[END_REF]. However, we believe that these analyses are not entirely complete and/or have certain drawbacks. For example, in [START_REF] Han | Probing Light Higgsinos in Natural SUSY from Monojet Signals at the LHC[END_REF] the 95% confidence level (CL) reach for the 14 TeV run of the LHC was calculated assuming a signal (S) to background (B) ratio below the 2% level, which is probably not quite feasible when taking into account that the actual systematic error should be above 3-5% even for quite optimistic analyses [START_REF] Khachatryan | Search for dark matter, extra dimensions, and unparticles in monojet events in proton-proton collisions at √ s = 8 TeV[END_REF][START_REF]Search for New Phenomena in Monojet plus Missing Transverse Momentum Final States using 10fb-1 of pp Collisions at sqrts=8 TeV with the ATLAS detector at the LHC[END_REF][START_REF]Sensitivity to WIMP Dark Matter in the Final States Containing Jets and Missing Transverse Momentum with the ATLAS Detector at 14 TeV LHC[END_REF]. In [START_REF] Han | Hunting Quasi-Degenerate Higgsinos[END_REF] the authors have performed their analysis at the parton level while, as shown in a preliminary analysis [START_REF] Brooijmans | Physics at TeV Colliders: New Physics Working Group Report[END_REF][START_REF] Belyaev | presentation at SUSY 2014 conference[END_REF], a fast detector simulation analysis leads to qualitatively different results and, therefore, is crucial. In [START_REF] Baer | Monojets and mono-photons from light higgsino pair production at LHC14[END_REF] the conclusion about the observability of the quasi-degenerate higgsino NSUSY scenario from the monojet search was negative. However in this study the authors did not attempt to optimise the E miss T which turns out to be important as we will show in this paper. In [START_REF] Schwaller | Compressed electroweakino spectra at the LHC[END_REF] the prospects were more optimistic even after including systematics uncertainties. However, also there no optimisation of the cuts was performed which considerably enhances the accessible mass range as we will show. One should also note the Ref. [START_REF] Baer | Monojet plus soft dilepton signal from light higgsino pair production at LHC14[END_REF] where authors suggested a new promising signature including a pair of soft leptons, and have demonstrated its potential power. However, in this study the important b b background was not considered as we discuss below, implying thus further background investigation for this signature which we do not consider at present. Finally one should mention Ref. [START_REF] Han | Accessing the core of naturalness, nearly degenerate higgsinos, at the LHC[END_REF], which studied similar to [START_REF] Baer | Monojet plus soft dilepton signal from light higgsino pair production at LHC14[END_REF] monojet plus soft lepton signature suggesting visibly harder cuts to suppress b b background. The respective higgsino mass reach from this study is quite limited. On the other hand the suggested b-jet veto will not quite work for the signature with isolated soft muons, so, we believe that one should estimate b b background more precisely even for the case of harder cuts suggested in [START_REF] Han | Accessing the core of naturalness, nearly degenerate higgsinos, at the LHC[END_REF].

Motivated by the above-mentioned previous studies, here we aim to perform a comprehensive and realistic analysis of the monojet potential to probe this NSUSY scenario. Our analysis is performed at the level of a fast detector simulation and the whole twodimensional NSUSY parameter space mentioned above, rather than selected benchmarks as were attempted previously, is explored therefore completely covering the region of our interest. We consider prudent systematic errors and optimise the kinematic cuts to keep the S/B ratio at a reasonable level. We then discuss the LHC potential to cover the NSUSY parameter space at 8 TeV and produced projections for the 13 TeV run of the CERN machine. By analysing also the exclusion potentiality for DM direct detection experiments, we aim to show that collider and DD experiments have a high degree of complementarity.

The paper is organised as follows. In Section 2 we describe the parameter space and mass spectrum of NSUSY while in Section 3 we discuss the DM properties of this scenario. Section 4 is dedicated to the analysis of the collider phenomenology of the compressed higgsino scenario while in Section 5 we show the complementarity of collider and the DD experiment. We conclude in Section 6.

Parameter space and Spectrum of NSUSY

In the bases ( B0 , W 0 , H0 d , H0 u ) and ( W 0 , H0 d ) the mass matrices of the neutralino and chargino sector of the MSSM are

M χ0 0 =      M 1 0 -M Z s ω c β M Z s ω s β 0 M 2 M Z c ω c β -M Z c ω s β -M Z s ω c β M Z c ω c β -µ M Z s ω s β -M Z c ω s β -µ 0      M χ± 1 = M 2 √ 2M W s β √ 2M W c β µ
(2.1) where M 1 and M 2 are the soft susy breaking mass parameter for B and W , µ is the Higgsino mass parameter, c ω and s ω are cos and sin of the Weinberg angle, tan β = v 2 /v 1 is the ratio of the vacuum expectation values of two Higgs doublets, s β , c β are sin β and cos β, respectively and m Z , m W are the masses of the SM gauge bosons Z 0 and W ±

As a first step we consider scenarios where the χ0 1 and χ± 1 have a high Higgsino component with all the SUSY partners of the SM fermions with masses in the multi TeV range and, to first obtain a qualitative understanding of the spectrum, we then expand the corresponding mass eigenvalues in the limit |µ| |M 1 |, |M 2 | obtaining:

m χ0 1,2 ∓ |µ| ∓ m 2 Z 2 (1 ± s 2β ) s 2 ω M 1 + c 2 ω M 2 (2.2)
m χ± 1 |µ| 1 + α(m Z ) π 2 + ln m 2 Z µ 2 -s 2β m 2 W M 2 (2.3)
where we have defined s 2β = sin(2β)sgn(µ) and α is the electromagnetic structure constant.

In the case of χ± 1 we have also included the electromagnetic corrections as this shifts the mass by about 0.5 % which is indeed important in those cases where the mass splitting between χ± 1 / χ0 2 and χ0 1 is of the order of a few GeV. We would like to note, that in our numerical results we have included the complete one-loop corrections in the calculation of the masses as this changes the absolute mass scale by several per-cent. However, the mass splittings discussed below are hardly affected as the additional corrections changes those at most by O(α/(4π) ln(m 2 W /m 2 Z ) which is below the per-mile level. For positive (negative) µ the mass eigenstate with negative (positive) CP-eigenvalue is the lightest one. The mass splittings are given by

∆m 0 = m χ0 2 -m χ0 1 m 2 Z s 2 ω M 1 + c 2 ω M 2 (2.4) ∆m ± = m χ± 1 -m χ0 1 ∆m 0 2 + |µ| α(m Z ) π 2 + ln m 2 Z µ 2 (2.5)
where we have neglected corrections of the order 1/ tan β and (µ/M 1,2 ) 2 . In order to analyse scenarios where the mass splitting between χ± 1 and χ0 1 varies from the (quasi) degenerate regime up to the regime with larger mass splittings, we have chosen the following range of parameter space: µ = (100, 300) GeV M 1 = (µ, µ + 600) GeV and (-µ, -µ -600) GeV tan β = (5 -50), (2.6) fixing the value of M 2 = 2 TeV, which has the effect of decoupling χ 0 4 and χ ± 2 that will not be considered anymore in the following, along with the rest of SUSY spectrum, which is assumed to be decoupled.

For |M 1 | µ one can obtain a simple approximation for χ0 1,2,3 masses [83] which we confirm by numerical evaluation, indicating that in this parameter region χ0

1 and χ0 3 are strongly mixed bino-higgsino states whereas χ0

2 is essentially a higgsino-like state with only a small bino component. µ, is clearly shown in this plots where the mass of the lightest neutralino, the DM candidate, is also presented.

Red: m χ 2 0 -m χ 1 0 Blue: m χ 1 + -m χ 1 0 M 1 =μ M 1 =μ+600 GeV
m χ 1 0 [GeV] Δ M[GeV] Red: m χ 2 0 -m χ 1 0 Blue: m χ 1 + -m χ 1 0 M 1 =-μ M 1 =-μ-600 GeV
m χ 1 0 [GeV] Δ M[GeV] Red: m χ 2 0 -m χ 1 0 Blue: m χ 1 + -m
The nature of the neutralinos and charginos, as well as the small mass splitting between them, has a strong impact on their decay modes. In the case of pure higgsinos, the three body decays are dominated by virtual vector bosons, and due to the small mass differences the decays into third generation fermions are suppressed. Note that in the scenario where M 1 is close to |µ|, the off-shell lightest Higgs boson, h 0 , can also give sizeable contributions [START_REF] Bartl | Large Higgs boson exchange contribution in three-body neutralino decays[END_REF][START_REF] Djouadi | Three body decays of SUSY particles[END_REF]. Also, one should note that in the case

|M 1 | µ, m χ0 3 -m χ0 1 2(m χ0 2 -m χ0 1 )

and χ0

3 decays to the lightest chargino with the 50% probability while sharing about 25% decay to each of χ0

1 and χ0 2 . Three body decays in the limit of small mass separation are discussed in [START_REF] Simone | Pseudo-Dirac Dark Matter Leaves a Trace[END_REF], where an effective theory study of the pseudo-Dirac DM scenario [START_REF] Nelson | The Minimal supersymmetric model without a mu term[END_REF][START_REF] Hsieh | Pseudo-Dirac bino dark matter[END_REF][START_REF] Belanger | Dark Matter with Dirac and Majorana Gaugino Masses[END_REF] such as the higgsino-like was performed. In this limit the decay width does not depend on the overall neutralino mass, but just on the mass difference

μ [GeV] M 1 [GeV] m χ 1 0 and m χ 1 + -m χ 1 0 [GeV]
0 μ [GeV] M 1 [GeV] m χ 1 0 and m χ 1 + -m χ 1 0 [GeV]
Γ( χ± 1 , χ0 2 → f f χ0 1 ) = C 4 128π 3 ∆m 5 Λ 4 (2.7)
where Λ m W,Z,h 0 is the mass of the leading mediator, while ∆m is either m χ0

2 -m χ0 1 or m χ± 1 -m χ0 1 .
As an example for off-shell Z exchange and decay into leptons, the coefficient C is

C 4 1 4 g 4 c 4 ω (s 2 ω -1/2) 2 + s 4 ω (2.8)
and a similar expression occurs for the case of the off-shell W -decay. The proper decay length is very sensitive to the value of ∆m, and values below the GeV lead to displaced vertices, or collider-stable χ± 1 and χ0 2 . Indeed, for the decay χ0 2 → f f χ0 1 with an off shell Z exchange, the proper decay length is given by

L = cτ 0.025 cm ∆m 1 GeV -5
(2.9) which implies that for ∆m 0.1 GeV, χ0

2 would be collider stable. Similarly, for ∆m 1 GeV one could look for displaced vertices of order 100 µm. Note that the measured decay length would depend on the boost factor of the decaying neutralino to be taken into account as discussed in details in Ref. [START_REF] Simone | Pseudo-Dirac Dark Matter Leaves a Trace[END_REF]. While collider stable χ0

2 will contribute to the E miss T , long lived or collider stable χ± 1 will provide a clear signature in the detector. A χ± 1 with a long enough lifetime can be detected in the tracking detectors by identifying decays that result in tracks with no associated hits in the outer region of the tracking system as recently analysed by ATLAS [START_REF] Aad | Search for charginos nearly mass degenerate with the lightest neutralino based on a disappearing-track signature in pp collisions at (s)=8 TeV with the ATLAS detector[END_REF] and CM [START_REF] Khachatryan | Search for disappearing tracks in proton-proton collisions at √ s = 8 TeV[END_REF] collaborations. Both collaborations have obtained similar results, concluding on sensitivity for charginos with a lifetime between 0.1 ns and 100 ns and covering chargino mass up to 500 GeV which significantly surpass the reach of the LEP experiments. For ∆M 0.25 GeV the chargino could be a collider-stable charged particle [START_REF] Simone | Leptogenic Supersymmetry[END_REF], and bounds on such a situation arising from the 8 TeV run of the LHC can be estimated to be m χ± 300 GeV [START_REF] Aad | Searches for heavy long-lived charged particles with the ATLAS detector in proton-proton collisions at √ s = 8 TeV[END_REF]. On the other hand, we have found that for ∆M 0.4 GeV there is no limit on m χ± 1 > 100 GeV from the above LHC searches. Therefore, our task is to analyse the potential of the monojet search to cover the NSUSY parameter space with ∆M 0.4 GeV.

Dark Matter direct and indirect detection in the NSUSY parameter space

The results from Planck [START_REF] Ade | Planck 2013 results. XVI. Cosmological parameters[END_REF][START_REF] Ade | Planck 2015 results. XIII. Cosmological parameters[END_REF] (see also WMAP [START_REF] Hinshaw | Nine-Year Wilkinson Microwave Anisotropy Probe[END_REF]) have further decreased the error on the already very precise measurement of the dark matter relic density, Ω Planck DM h 2 = 0.1184 ± 0.0012.

As we assume R-parity to hold, the LSP will be stable and will contribute to this relic density. In the scenarios under consideration, the LSP is the lightest neutralino, χ0

1 , which is dominantly higgsino-like with a variable bino component. It is well known that for χ0

1 of mass O(100 GeV), a higgsino-like LSP provides DM relic density below the level observed by Planck. This happens because of the annihilation and co-annihilation rate of LSP and NLSP particles in the early Universe begin too large. On the other hand, for bino-like neutralinos the annihilation is suppressed, resulting in DM over-abundance. In the mass range we study, mixed bino-higgsino LSPs can therefore lead to the correct relic abundance. For µ M 1 the LSP is mainly higgsino-like and the value of Ω DM h 2 Ω Planck DM h 2 , however this nonetheless at least solves the typical problem of the over-closure of the universe for neutralinos of this mass range in generic SUSY parameter space. In this case we then assume that the remaining relic abundance is accounted for by other means, for example, it could come from multi-TeV moduli field where the higgsino LSP is non-thermally produced (see e.g. Ref. [START_REF] Allahverdi | Non-thermal Higgsino Dark Matter: Cosmological Motivations and Implications for a 125 GeV Higgs[END_REF]) or from mixed axion-higgsino DM (e.g. Ref. [START_REF] Baer | Natural Supersymmetry: LHC, dark matter and ILC searches[END_REF]).

In order to assess the compatibility of the scenarios under our investigation with existing experimental limits, we have evaluated Ω DM h 2 , the spin-independent annihilation cross section (σ SI ) and the respective DD rates using micrOMEGAs 2.4.1 [START_REF] Belanger | MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model[END_REF][START_REF] Belanger | Indirect search for dark matter with micrOMEGAs2.4[END_REF]. In Fig. 3 we show the results for Ω DM h 2 as a function of µ for tan β = 5, 50 and positive or negative values of M 1 . Note that the uncertainty on Ω DM h 2 is not shown, the full one-loop corrections are not yet available, but we expect that these are not too large and will not qualitatively change our conclusions. From these plots we see, as expected, that in general Ω DM h 2 lies below Ω Planck DM h 2 , and decreases as χ0 1 becomes increasingly higgsino-like. This is because of the mass splitting between χ± 1 and χ0 1 , which becomes larger as M 1 decreases. This suppresses the coannihilation channels which otherwise lead to an efficient reduction of the relic density. Below the W W threshold, the annihilation of the mixed gaugino-higgsino neutralinos therefore occurs via the Higgs and Z-bosons. The small bottom Yukawa coupling and the suppressed coupling to the Z implies that this mechanism is not efficient, apart from at the Z and h resonance. Therefore the spike in the relic density µ ∼ 100 GeV can be explained by the fact that this is just below the W W threshold. At lower values of µ ∼ 70, 90 GeV one moreover observes two dips corresponding to the Z-boson and Higgs funnels. This is most pronounced for positive values of M 1 and lower values of tan β where the mass splitting is larger.

In Fig. 4 we further show the spin-independent annihilation cross section for DD, again for positive and negative M 1 as in Fig. 3, where instead of σ SI we plot the rescaled quantity R Ω σ SI (pb), where the scaling factor R Ω = Ω DM /Ω Planck DM allows easy comparison with the most recent limits (also reported in these plots) from LUX [START_REF] Akerib | First results from the LUX dark matter experiment at the Sanford Underground Research Facility[END_REF], as well as the projected limits from XENON1T after 2 years live-time and 1 ton fiducial mass (see e.g. Ref. [START_REF] Aprile | The XENON1T Dark Matter Search Experiment[END_REF]), which in general assume the relic density to be the value measured by Planck. Fig. 4 illustrates that the region with a low LSP masses and a higher mass splitting between χ± 1 and χ0 1 (see Fig. 2), which is also easiest to see at colliders, is in fact excluded as the σ SI for the DD experiments is too high. In the following we will further highlight the interesting complementarity between the reach of the collider searches and the DD searches, particularly interesting for low DM masses. In Fig. 5 we also show the mass splitting between χ± 1 and χ0 1 in the µ-M 1 plane of Fig. 2, along with the region excluded by LUX and the projected exlusion regions from XENON1T. This emphasises that the region where µ and M 1 are very close each other is already excluded by LUX, which in turns puts an upper bound on the splitting between χ± 1 or χ0 2 and χ0 1 . Note that for positive values of M 1 , the splitting is larger, as is the mixing between bino and higgsino component of the DM. As larger mixing leads to larger couplings to Higgs bosons, DD is more sensitive to the case M 1 > 0.

A related question is whether these scenarios could be excluded by indirect detection (ID) experiments, i.e. the detection of energetic e ± , γ, p or p, which may be created by the pair annihilation of weakly interacting massive particles (WIMPS). It turns out that the strongest bounds on neutralinos coming from such experiments are set by gamma ray telescopes: both the Fermi-LAT gamma-ray space telescope [START_REF] Ackermann | Constraining Dark Matter Models from a Combined Analysis of Milky Way Satellites with the Fermi Large Area Telescope[END_REF] as well as ground based telescopes. Fermi-LAT is sensitive to gamma rays particularly in the low mass range up to O(100 GeV). It is therefore particularly sensitive to lighter mixed bino-higgsino neutralinos, but the bounds are not competitive with those coming from DD. Here we are interested in light higgsinos with possibly some bino component, for which the relic density is in general much below the value measured at Planck/WMAP (see Fig. 3). In this case, all bounds from ID must be scaled by the square of the ratio of the predicted relic density to the experimental value. The rescaling appears as the pair annihilation cross section depends on the square of the local WIMP abundance. In [START_REF] Baer | Direct and indirect detection of higgsino-like WIMPs: concluding the story of electroweak naturalness[END_REF][START_REF] Arrenberg | Working Group Report: Dark Matter Complementarity[END_REF] it was shown that the Fermi-LAT limits derived for WIMP annihilations into W W (of which a large component of the total annihilation cross section should be comprised) are not yet sensitive to higgsino LSPs in the 100-350 GeV range mainly due to the predicted under-abundance of neutralinos in these scenarios. Note that WIMPS having masses in the range 200 GeV to a few TeV will be probed by the future CTA array [START_REF] Actis | Design concepts for the Cherenkov Telescope Array CTA: An advanced facility for ground-based high-energy gamma-ray astronomy[END_REF][START_REF] Doro | Dark Matter and Fundamental Physics with the Cherenkov Telescope Array[END_REF], i.e. any higgsino or wino-like LSP for which the relic abundance is within an order of magnitude of Planck would be seen. In [START_REF] Cahill-Rowley | Complementarity of dark matter searches in the phenomenological MSSM[END_REF] a comprehensive scan of the pMSSM was carried out, and DM limits (both present and projected) coming from the LHC, DD and ID experiments were studied (note that the collider study was not dedicated to the region of interest of the present paper; our optimisation of the kinematical cuts goes beyond previous analyses). Ref. [START_REF] Cahill-Rowley | Complementarity of dark matter searches in the phenomenological MSSM[END_REF] found an impressive degree of complementarity: in the range a few hundred GeV to 1 TeV DD poses stronger bounds, and for higher masses ID is more sensitive. Therefore we conclude that for the mass range studied in this paper, the LHC and DD limits are the most relevant.

LHC potential to probe NSUSY

In the previous section we have discussed the current and future sensitivity of underground experiments to NSUSY scenario. In this section we explore the LHC potential to probe NSUSY and demonstrate that it plays a crucial complementary role.

The LHC's most sensitive searches for quasi-degenerate χ± i and χ0 i scenario are monojet signatures, i.e. the production of a pair of electroweakinos through the s-channel exchange of a SM EW gauge boson, Z, γ or W ± , accompanied by hard QCD initial state radiation (ISR) via the process

pp → χ a χ b j χ a,b = χ0 1,2,3 , χ± 1 . (4.1)
The corresponding Feynman diagrams can be found in Fig. 6.

q(q ′ ) q Z(W ± ) q g χ0 1,2,3 χ0 1,2,3 (χ ± 1 ) q g q q(q ′ ) Z(W ± ) χ0 1,2,3 χ0 1,2,3 (χ ± 1 ) q(q ′ ) q Z(W ± ) g q χ0 1,2,3 χ0 1,2,3 (χ ± 1 ) q q Z/γ q g χ- 1 χ+ 1 q g q q Z/γ χ- 1 χ+ 1 q q Z/γ g q χ- 1 χ + 1 Figure 6.
Representative diagrams for pair neutralino-chargino production in association with quark/gluon leading to monojet signature.

As we discuss below, the main problem for the signal search in this channel is the large background, the dominant contribution coming from the EW production processes Z+jets and W +jets. We begin our analysis of the monojet signal from the recasting of the LHC at √ s = 8 TeV data and the respective experimental results, and then analyse the prospects of the 13 TeV LHC run with both standard and high luminosity (HL) options.

Both ATLAS [START_REF]Search for New Phenomena in Monojet plus Missing Transverse Momentum Final States using 10fb-1 of pp Collisions at sqrts=8 TeV with the ATLAS detector at the LHC[END_REF] and CMS [START_REF]Search for new physics in monojet events in pp collisions at sqrt(s)= 8 TeV[END_REF] have performed studies of monojet signatures at the LHC Run1, which have been interpreted in the context of an EFT approach. However this approach cannot be used for the NSUSY scenario because of the Z, W, γ mediating interactions as indicated in Fig. 6. Hence, in the case of NSUSY,

DM with EW mediators / ⇒ L ef f = 1 Λ 2 (qΓq) (χ χ) (4.2)
where Γ is some Lorentz structure and χ is the DM particle. Therefore it is necessary to recast the searches in terms of this NSUSY scenario and not to use limits from Λ-DM mass plane.

Analysis Setup

In this section we describe the different aspects of our simulation of both the signal and the most important backgrounds, implementing the important steps of hadronization and fast detector simulation. We performed a parton-level simulation using MadGraph v1.5.11 [START_REF] Alwall | MadGraph 5 : Going Beyond[END_REF] with the MSSM model available on the FeynRules web page [START_REF] Christensen | FeynRules -Feynman rules made easy[END_REF] implemented in UFO format [START_REF] Degrande | UFO -The Universal FeynRules Output[END_REF] 1 , and cross-checked results against CalcHEP [START_REF] Belyaev | CalcHEP 3.4 for collider physics within and beyond the Standard Model[END_REF] with the MSSM model from the HEPMDB website 2 .

At the level of matrix-element we have generated the production of a pair of electroweakinos via the s-channel exchange of a SM EW gauge bosons, accompanied by hard QCD initial state radiation. Parton level SM background simulations have been also cross checked between two packages. Our choice of PDF sets is CTEQ6L1 [START_REF] Pumplin | New generation of parton distributions with uncertainties from global QCD analysis[END_REF] and we used the MadGraph dynamical choice of renormalization scale which is equal to the geometric mean of Mass 2 + P 2 T for the final state particles. Parton showering, hadronisation and decay of the unstable particles were simulated using PYTHIA v6.4 [START_REF] Sjostrand | PYTHIA 6.4 Physics and Manual[END_REF] while detector effects have been simulated with Delphes3 [START_REF] De Favereau | DELPHES 3, A modular framework for fast simulation of a generic collider experiment[END_REF] employing a suitable CMS card. Finally, the background processes yields, which include Z + j, W + j, t t, QCD and single top processes, have been taken from the experimental results for 8 TeV analysis and has been simulated using MadGraph+PTYHIA+DELPHES chain for the 13 TeV analysis.

For the signal and background events at the parton level we have applied a cut on the jet transverse momentum of p j T > 90 GeV at the generation level. This is a subtle point. One might argue that, since dealing with processes involving QCD radiation, it is necessary to apply a merging procedure between the hard jet generated via the parton-level matrix element and the soft jets generated by the showering algorithm, merging therefore the 0-jet and 1-jet samples. However, since a hard jet is selected at the analysis level (e.g. p j T > 110 GeV for the 8 TeV LHC analysis and higher for the 13 TeV case, with a final selection requirement of a high E miss T , somewhat correlated with the jet p T [START_REF] Brooijmans | Physics at TeV Colliders: New Physics Working Group Report[END_REF]) we found that this matching was unnecessary, and generated just the one jet sample. Moreover, avoiding matching in this case allowed us to make our analysis much more effective and obtain enough statistics in the high p j T region. To illustrate the validity of this procedure, in Fig. 7 we present the p T distribution of the leading jet for the case of χ0 Shown in blue is the leading jet p T for the case of a 0+1 matched jet sample, while the case of the 1-jet events being unmatched is shown in black. One can indeed see that in case of a high p T cut on the leading jet at the analysis level, the matched and unmatched distributions of the leading jet are very similar, indicating that contribution from the 0-jet matched sample is negligible, and that the p T of the leading jet is dominated by the one-jet sample in the high p T region.

LHC Run1: the reach of monojet searches

We start our analysis with the exploration of the 8 TeV LHC potential to probe NSUSY and the recasting of the respective experimental results. For this purpose we have chosen CMS monojet analysis [START_REF]Search for new physics in monojet events in pp collisions at sqrt(s)= 8 TeV[END_REF] which has been done for the data recorded at √ s = 8 TeV with 19.5 fb -1 integrated luminosity3 . We have applied the following trigger selection followed by the cut-flow according to the CMS analysis:

• Two triggers which require E miss T > 120 GeV or E miss T > 105 GeV and a jet with p T > 80 GeV and within |η| < 2.6.

• The analysis then requires that the jet with the highest transverse momenta has p T > 110 GeV and |η| < 2.4.

• Events with more than two jet with p T > 30 GeV and |η| < 4.5 are discarded together with events where ∆φ(j 1 , j 2 ) < 2.5, where j 1 and j 2 are the leading and sub-leading jets, to reduce QCD background.

• The W production background was suppressed by applying a veto on events with one electron or muon satisfying p T > 10 GeV cut and events with one tau jet with p T > 20 GeV and |η| < 2.3.

• Finally the analysis was performed in 7 regions with an increasing requirement of E miss T : E miss T > 250, 300, 350, 400, 450, 500 and 550 GeV.

We have then derived the signal significance for each signal region from the number of signal events (S), after having imposed the cuts above, and the number of background events (B), using the following expression as in the CMS analysis:

α 2( √ S + B - √ B) , (4.3)
which is similar to the more common S/ √ S + B for S B, but is more robust with respect to downward fluctuations [START_REF] Bityukov | On observability of signal over background[END_REF]. In Fig. 8 we show contours of iso-significance α (solid gray) in the µ-M 1 plane for the lowest (left frame) and highest (right frame) requirements on E miss T . The LUX exclusion is further shown by the green shaded area as well as the LEP2 limit on charginos by the red shaded area. In addition, the S/B ratio is shown by a blue-dashed line. By only inspecting significance contours, one could think that the Run1 LHC data could slightly extend the LEP2 limit using the low E miss T signal region: the area on the left of the lines of α=2 would be ruled out at 95% confidence level (CL). However, one then observes that the S/B never goes above 2% in the region allowed by LEP2, while the actual systematic uncertainties of this analysis are of the order of 5-10% [START_REF]Search for New Phenomena in Monojet plus Missing Transverse Momentum Final States using 10fb-1 of pp Collisions at sqrts=8 TeV with the ATLAS detector at the LHC[END_REF][START_REF]Search for new physics in monojet events in pp collisions at sqrt(s)= 8 TeV[END_REF]. This means that the low S/B ratio for the 8 TeV LHC is the main obstacle to going beyond the LEP2 limits. We can also see from Fig. 8 (right) that the higher E miss T cut increases the S/B ratio, eventually at the expense of the signal, meaning that the significance drops below α=2 level in the parameter space allowed by LEP2. Therefore we conclude that Run1 LHC does not have the potential to test the NSUSY scenario beyond the LEP2 limit due to systematic uncertainties.

13 TeV LHC potential and complementarity to underground experiments

As we have seen in the previous section, the LHC Run1 is not sensitive to the NSUSY parameter space we consider due to the low statistical significance and the fact that the low S/B ratio remains below the systematic errors.

In this section we study the LHC Run2 case and show that the higher collider luminosity and energy allow us to choose kinematical cuts, bringing the S/B ratio to a desirable level while keeping the statistical significance at a high enough level in order to establish sensitivity to the NSUSY parameter space. Being a very challenging scenario, we make projections up to HL configuration of the LHC machine for different assumptions regarding the E miss T cuts and the control over systematic uncertainties. In order to ensure that the S/B ratio is under control, we compare the relative size and shape difference of the signal versus the dominant irreducible background Z +jet → ν ν+jet (Zj). The relevant parton level distribution are shown as a function of the jet p T in Fig. 9. One can see that even for µ = 93 GeV corresponding to m χ 0 1 100 GeV, the background is about 3 orders of magnitude higher than the signal for a low p j T > cut. An important feature of the signal versus background is that the shape of the background distribution is quite different from the signal: the background falls more rapidly with p j T , and the difference in the slope with respect to the signal is bigger for higher neutralino masses. The different in slope is mainly due to the mass difference between the neutralino, from the signal, and the neutrino, from the background. One should also notice that the difference between the shapes of the signal and background p j T distributions vanishes for very large values of p j T m χ 0 1 , as one would expect. This provides us with our main strategy to optimise the LHC sensitivity to the NSUSY parameter space, that is to find the maximal value for E miss T cut while maintaining the statistical significance at high enough level. On inspecting this figure one can already see that in order to achieve S/B ratio at about 5% level the value E miss T cut should be around of 1 TeV. The signal process analysed is the same as that in Eq. (4.1), with the only difference being that the initial requirement on the QCD ISR has been increased to 150 GeV, motivated by our preliminary study in Ref. [START_REF] Brooijmans | Physics at TeV Colliders: New Physics Working Group Report[END_REF].

Along with the signal processes, we have simulated the two main backgrounds 1. Cutflow for the two main SM background and two choices of signal for the 13 TeV LHC with 100 fb -1 of integrated luminosity. The initial number of events corresponds to p j T > 150 GeV cut.

pp → Zj → ν l νl j pp → W j → l + ν l + c.c. ( 4
Zj, Z → νν W j, W → lν µ = 100 GeV µ = 200 GeV M 1 =700 GeV M 1 =800 GeV Initial #
• We apply a veto on events with more than two jets with p T > 30 GeV and |η| <4.5

• We require ∆φ(j 1 , j 2 ) <2.5

• We apply a veto on electrons and muons with p T > 10 GeV.

• We apply a veto on taus with p T > 20 GeV and |η| < 2.3.

We have then defined signal regions with increasing cuts on E miss T ; an example cut flow is provided in Table 1. We have not generated t t, QCD and single top background, though we have applied to our simulated samples cuts that reduce these background to a negligible level with respect to W j and Zj, see details in Refs. [START_REF]Search for New Phenomena in Monojet plus Missing Transverse Momentum Final States using 10fb-1 of pp Collisions at sqrts=8 TeV with the ATLAS detector at the LHC[END_REF][START_REF]Search for new physics in monojet events in pp collisions at sqrt(s)= 8 TeV[END_REF].

As observed for the case of 8 TeV LHC, a strong tension arises in attempting to simultaneous maximise the S/B ratio and the statistical significance. We demonstrate this in Fig. 10, where we plot S/B and α as a function on the final selection cut on the E miss T for both the case of χ0

1 ∼ 100 GeV (left panel) and ∼ 200 GeV (right panel). This figure clearly indicates that in order to achieve a high enough S/B ratio in keeping with the expected level of systematic uncertainties, a hard cut on E miss T ought to be applied, which at the same time pulls the significance down, due to the reduction in the number of signal events. Given this tension it is therefore important to optimize the E miss T cut to provide a high enough S/B ratio and to keep α >2 (5) in order to obtain an exclusion (discovery). The NSUSY parameter space can be conveniently described and presented in m χ0 1 -∆M plane, where for ∆M we have chosen the mass difference between the lightest chargino and the DM particle. In this plane, for a given value of integrated luminosity, the optimal E miss T cut can be chosen by the point where S/B and α cross or are as close to each other as possible. This is related to the fact that the iso-significance contours are shifted to the left in the m χ0

1 -∆M plane with the increase in the E miss T cut due to the decrease of signal statistics, while iso-S/B contours are shifted to the right at the same time due to the increase of S/B ratio. Therefore the case when the respective iso-contours cross/are close to each other, would provide the maximal exclusion or discovery area in the m χ0 1 -∆M plane. We illustrate this in Fig. 11 which presents S/B and significance isocontours in the m χ0

1 -∆M plane for two different cuts on E miss T , 850 and 900 GeV. One can see that indeed for E miss T > 850 GeV, the exclusion area is below S/B = 3% (blue dashed) contour, while for E miss T > 900 GeV, the area below α = 2 (red solid) contour is excluded. Since for the first case the exclusion area is bigger, the E miss T > 850 GeV is better choice for the optimal cut. We have found that a cut around 600 (850) GeV for 100 fb -1 (3000 fb -1 ) provides α 2 and S/B 0.03 iso-contours optimally close to each other, which maximises the reach of the 13 TeV LHC for the NSUSY parameter space. The proximity of α 2 and S/B 0.05 iso-contours eventually requires higher E miss T cut which is found to be around 950 GeV, and as a result leads to a poorer 13 TeV LHC reach as we discuss below. While presenting results in the m χ0 1 -∆M plane, we separate the cases M 1 > 0 and M 1 < 0, which differ due to different bino-higgsino components and mass splittings as discussed earlier, which eventually becomes less and less relevant as we approach the low ∆M region.

In Fig. 12 and Fig. 13 we show results for M 1 > µ and M 1 < -µ (note that we have chosen µ > 0). The left panels contain the 2σ exclusion LHC reach while the right panels contain the 5σ discovery LHC potential. Both the cases of requiring 3% and 5% for the S/B ratio are shown, the latter just for the High Luminosity (HL)-LHC scenario. In the same plot we present also LUX sensitivity and the projected sensitivity of XENON1T to the in the m χ0 1 -∆M plane. One can see that at Run2, with ∼ 100 fb -1 of integrated luminosity collected, the LHC will be able to exclude up to ∼150 GeV χ0 1 with ∆M below 5 GeV, if the systematic error can be kept at the 3% level. The LHC will therefore surpass the LEP sensitivity for this scenario which reached the limit χ± 1 > 103 GeV for low ∆M values as shown in Fig. 12 and Fig. 13. By the end of the HL-LHC run, up to ∼ 250 (200) GeV LSP could be excluded for low mass splitting with a S/B> 3% (5%). It is quite remarkable that the LHC has the maximal sensitivity in the low ∆M region. This reach is nicely complemented by the LUX results and the projected exclusions for XENON1T, which cover the region with higher mass splitting. Such complementarity would allow LHC and DD experiments to completely exclude NSUSY scenario with DM mass up to ∼ 250 GeV. Discovery prospects for this scenario are shown in the right panel of Fig. 12 which demonstrates that ∼ 180 (110) GeV χ 0 1 can be discovered for S/B > 5% (3%) at the end of the HL-LHC run, while with 100 fb -1 the LHC will not have discovery sensitivity. At the end of this section we would like to discuss a potentially promising new signature involving a monojet plus soft di-leptons which was studied in terms of the NSUSY parameter space in Ref. [START_REF] Baer | Monojet plus soft dilepton signal from light higgsino pair production at LHC14[END_REF]. In short, in this paper the signal from the second neutralino decaying leptonically into the lightest neutralino has been studied. It was suggested that one could trigger on soft leptons and use the upper cut on di-lepton invariant mass below 10 GeV to suppress background and extract the signal from compressed χ0

1 -χ0 2 production. After the suggested cuts, the background can be reduced down to about the 6 fb level bringing it below the signal, which would allow one to claim a discovery for chosen benchmarks with only 100 fb -1 of integrated luminosity. One should note however, that in spite of the comprehensive set of backgrounds, jb b was not considered in this paper. After application of the cuts from [START_REF] Baer | Monojet plus soft dilepton signal from light higgsino pair production at LHC14[END_REF] we have found that jb b background is in fact dominant for the monojet plus soft di-lepton signal and is about two orders of magnitude above the backgrounds taken into account in the aforementioned study. An estimation of jb b background is not trivial, since the cross section of the jb b process is very high, about 100 nb for soft initial cuts on jet p T , while the efficiency of the selection cuts is very low, about 10 -5 . Therefore, a reasonable estimate of this background requires either the simulation of at least 10 6 events or the extraction of this background from future experimental data. Our preliminary results indicate that one should develop a dedicated strategy to suppress this jb b background in order to make the monojet plus soft di-lepton signal a viable tool ( 3 % ) LHC13 3 ab -1 (5%) Figure 13. Exclusion (left) and discovery (right) contour lines for the 13 TeV LHC at the end of the LHC Run2 (light red region) and of the HL-LHC (light blue region) assuming S/B>3%. For the latter case also the case S/B>5% is shown. The region excluded by LUX and the projected exclusion by XENON1T are also shown, together with the LEP limit on the χ± 1 mass. M 1 < -µ is considered here.

for the exploration of the NSUSY parameter space.

Conclusions

In this paper we have explored the complementary potential of the Large Hadron Collider and underground experiments to probe Dark Matter (DM) in the Natural Supersymmetry (NSUSY) scenario. This study, which combines searches from different kinds of experiments, has to be done in the context of a specific model, as (model-independent) Effective Theory (EFT) approaches are very limited in scope, see e.g. the discussion in Refs. [START_REF] Buchmueller | Beyond Effective Field Theory for Dark Matter Searches at the LHC[END_REF][START_REF] Busoni | On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC[END_REF]. In particular the EFT approach is not applicable for well motivated NSUSY scenario, which we study here, where DM has direct couplings to Standard Model electroweak (EW) gauge bosons and the Higgs.

Current limits on simple SUSY scenarios are at the TeV range, in clear tension with naturalness arguments and hence with the motivation for introducing SUSY in the first place. A possible explanation for this situation is that the manifestation of SUSY is not as simple as one expects, but there is more complexity in the structure of SUSY at highenergies. Notwithstanding, one would still expect that the particles more directly related to the tuning of the EW scale remain light in the spectrum. This leads to a generic expectation that DM in NSUSY should have a sizeable Higgsino component.

While being theoretically attractive this scenario also represents a clear example of how colliders and underground experiments can complement each other. Indeed, while underground experiments have a larger mass sensitivity than colliders, being able to probe the multi-TeV region, colliders can cover parameter space hidden from DM direct detection (DD) experiments. Specifically, the increase of the DM higgsino component makes NSUSY parameter space increasingly difficult to probe in DD experiments. In this region higgsinolike DM is quasi-degenerate with two other particles, the second neutralino and the lightest chargino, and the increase of higgsino component is correlated with the decrease of this mass splitting, ∆M. At the same time, the LHC sensitivity increases with the increase of the higgsino component of DM and reaches its maximum for very low values of ∆M. We have conveniently described the NSUSY parameter space in the µ-M 1 region and have translated it into sensitivity of the LHC and DM DD experiments in the m χ0 1 -∆M plane, where for ∆M we have chosen the mass difference between the lightest chargino and DM. We have studied the current and the future reach of underground experiments in this region, as well as given predictions for the relic abundance in NSUSY.

We would like to point out that we have presented combined LHC and DM DD experiments results for the whole m χ0 1 -∆M NSUSY space, rather than for chosen benchmarks. Moreover, we have optimised the final E miss T selection cut to keep the S/B ratio at a high enough level, to deal with the systematic errors on the background, which aer one of the main problem for the exploration of the NSUSY parameter space. As a result, we have found that the 8 TeV LHC unfortunately does not allow to test the NSUSY parameter space beyond the LEP2 limits, while the 13 TeV LHC at 3 ab -1 has the potential to significantly surpass the LEP2 reach and to cover DM masses up to about 250 GeV for ∆M<5 GeV, which can neither be covered by LUX nor by XENON1T DM DD experiments. At the same time the XENON1T experiment will be able to complementarily cover the ∆M>5 GeV parameter space up to large values of the DM mass, well beyond the LHC reach via monojet analyses.

Figure 1 . 1 -m χ0 1

 111 Figure 1. χ± 1 -χ0 1 and χ0 2 -χ0 1 mass splitting values as a function of m χ 0 1 (upper row) and M 1 (lower row) for the case M 1 > 0 (left) and M 1 < 0 (right).

Figure 2 .

 2 Figure 2. Mass splitting between χ± 1 and χ0 1 (blue solid line) and mass of χ0 1 (gray dashed line) in the µ-M 1 plane.

Figure 3 .

 3 Figure 3. The predicted value of the dark matter relic density Ω DM h 2 is shown as a function of µ for tan β = 5, 50 and positive (left) or negative (right) values of M 1 as indicated. The relic density measured by the Planck satellite, Ω Planck DM h 2 , is also shown for comparison, and the region excluded due to an overabundance of DM is indicated in grey.

Figure 4 .

 4 Figure 4. The logarithm of the predicted value of the spin-independent annihilation cross section for DD σSI = R Ω σ SI /(10 -8 pb), rescaled by R Ω /(10 -8 pb) where R Ω = Ω DM /Ω Planck DM , is shown as a function of µ for tan β = 5, 50 and positive (left) or negative (right) values of M 1 as indicated. The excluded limit from LUX (green), as well as the projected exclusion from XENON1T (cyan) are also shown for comparison, where the solid and dashed lines represent the exclusions for |M 1 | = µ and |M 1 | = µ+600 GeV respectively.

Figure 5 .

 5 Figure 5. The mass splitting between χ± 1 / χ0 1 (blue solid line) and mass of χ0 1 (black dashed line) is shown in the µ-M 1 plane, along with the region excluded by LUX results (green).

Figure 7 . 1 χ0 2

 712 Figure 7. Leading jet p T cross section distributions for the case of the 0+1 jet sample (blue) and 1 jet sample (black) for the χ0 1 χ02 production at the 8 TeV LHC.

Figure 8 .

 8 Figure 8. Contours for the 2( √ S + B -√ B) (gray) and S/B (blue-dashed)in the plane (µ,M 1 ) for the signal regions 1 and 7 as defined in [105]. LUX and LEP exclusions are shown in green and red.

. 4 )

 4 with l = e, µ, τ .Inspired by the 8 TeV CMS monojet search we have then applied the following cut-flow• We require a leading jet with p T > 200 GeV and |η| < 2.4

Figure 9 .

 9 Figure 9. Signal (dotted blue and dashed red) and Zj background (solid black) parton-level p j T distributions for the 13 TeV LHC for the NSUSY scenario. Left: p j T distributions for 100 fb -1 integrated luminosity. Right: normalised signal and Zj background distributions.

m χ 1 0m χ 1 0 1 Figure 10 .χ0 1 ∼

 111101 Figure 10. S/B and α as a function on the final selection cut on the E miss T for χ0 1 ∼ 100 GeV (left panel) and ∼ 200 GeV (right panel).

Figure 11 .

 11 Figure 11. S/B (blue) and α (red) isocontours for two choices of E miss T cut: 850 GeV (dashed) and 900 GeV (solid) in the m χ 0 1 -∆M plane.

Figure 12 .

 12 Figure 12. Exclusion (left) and discovery (right) contour lines for the 13 TeV LHC at the end of the LHC Run2 (light red region) and of the HL-LHC (light blue region) assuming S/B>3%. For the latter case also the case S/B>5% is shown. The region excluded by LUX and the projected exclusion by XENON1T are also shown, together with the LEP limit on the χ± 1 mass. M 1 > µ is considered here.

Table

		of events	3.15•10 6	1.25•10 7	3.63•10 5	6.45•10 3
	p j T > 200 GeV |η j | < 2.4 Jet veto	1.05•10 6 8.7•10 5	4.11•10 6 3.13•10 6	1.73•10 5 1.33•10 5	3528 2691
	∆φ(j 1 , j 2 ) < 2.5 Veto e ± , µ ± , τ ± E miss T > 200 GeV	7.2•10 5 7.2•10 5 6.4•10 5	2.3 •10 6 6.8 •10 5 4.3•10 5	1.10•10 5 1.08•10 5 9846	2320 2301 2188
	E miss T	> 600 GeV	4353	1002	171	93
	E miss T	> 700 GeV	1703	250	80	47
	E miss T	> 800 GeV	694	0	37	22

http://feynrules.irmp.ucl.ac.be/wiki/MSSM

http://hepmdb.soton.ac.uk/hepmdb:0611.0028

The respective ATLAS analysis, which leads to similar results, can be found in[START_REF] Aad | Search for new phenomena in final states with an energetic jet and large missing transverse momentum in pp collisions at √ s = 8 TeV with the ATLAS detector[END_REF].

Acknowledgements

The work of VS and AB is supported by the Science Technology and Facilities Council (STFC) under grant number ST/L000504/1 and ST/L000296/1 respectively. WP is supported by the Bundesministerium für Bildung und Forschung (BMBF) under contract no. 05H12WWE. AKMB thanks Andreas Goudelis for helpful discussions. AB and DB thank Marc Thomas for useful discussions and help with solving of jet matching problems.


Komi Soho
email: akomi.soho@univ-lorraine.fr

Farid Abed-Meraim
email: bfarid.abed-meraim@ensam.eu

Xavier Lemoine
email: cxavier.lemoine@ensam.eu

Hamid Zahrouni
email: dhamid.zahrouni@univ-lorraine.fr

Multiscale finite element simulation of forming processes based on crystal plasticity

Keywords:

For the numerical simulation of sheet metal forming processes, the commercial finite element software packages are among the most commonly used. However, these software packages have some limitations; in particular, they essentially contain phenomenological constitutive models and thus do not allow accounting for the physical mechanisms of plasticity that take place at finer scales as well as the associated microstructure evolution. In this context, we propose to couple the Abaqus finite element code with micromechanical simulations based on crystal plasticity and a selfconsistent scale-transition scheme.

This coupling strategy will be applied to the simulation of rolling processes, at different reduction rates, in order to estimate the evolution of the mechanical properties. By following some appropriately selected strain paths (i.e., strain lines) along the rolling process, one can also predict the texture evolution of the material as well as other parameters related to its microstructure. Our numerical results are compared with experimental data in the case of ferritic steels produced by ArcelorMittal.

Introduction

Although commonly used for the simulation of sheet metal forming processes, the traditional commercial finite element (FE) codes have several limitations. These restrictions are mainly related to the built-in constitutive models, which are essentially phenomenological thus not accounting for the physical mechanisms of plasticity and the underlying microstructure evolution.

To overcome some of these limitations, an alternative would be to couple these FE simulation codes with micromechanics-based behavior models. The latter are generally based on crystal plasticity, for the single crystal constitutive equations, and a scale-transition scheme (e.g., the selfconsistent model), for the overall response of the polycrystalline aggregate.

In the literature, several attempts have been made in order to implement constitutive models based on crystal plasticity (strain-rate dependent or strain-rate independent) into each integration point of FE calculations.

For strain-rate dependent constitutive models, we can mention, among others, the recent works by Segurado et al. [START_REF] Segurado | Multiscale modeling of plasticity based on embedding the viscoplastic self-consistent formulation in implicit finite elements[END_REF], who embedded the visco-plastic self-consistent formulation in implicit finite elements and by Zhang et al. [START_REF] Zhang | An effective semi-implicit integration scheme for rate dependent crystal plasticity using explicit finite element codes[END_REF], who proposed an effective semi-implicit integration scheme for rate dependent crystal plasticity using explicit FE codes. For strain-rate independent behavior models, several contributions can be found in the literature, among which the works by Scacciatella [START_REF] Scacciatella | Intégration d"un modèle de comportement élastoplastique à transition d'échelle dans un code de calcul de structures par éléments finis[END_REF], in the framework of small-strain elasto-plasticity, and by Zattarin [START_REF] Zattarin | Etude de l"intégration d'un modèle polycristallin dans un code d'élements finis en élastoplasticité[END_REF], in the large-strain framework.

However, the main drawback of such direct coupling approaches is the computation time that is relatively high. In fact, the CPU time in a polycrystalline model is generally proportional to the number of constituent grains (it also depends on the constitutive law adopted for the crystallographic slip), and this model is called at each integration point. Accordingly, the CPU time required for FE calculations based on a polycrystalline model is extremely high as compared to that associated with a phenomenological model.

For the above-mentioned reasons, another strategy is proposed in the current work, in which the polycrystalline model and the FE simulation codes Abaqus and LAM3 will be coupled in an indirect way. This technique will be applied to simulate rolling processes, at different reduction rates, and to estimate the evolution of the mechanical properties. By following some appropriately selected strain paths (i.e., strain lines) along the rolling process, we are also able to predict the texture evolution of the material as well as other parameters related to its microstructure.

Constitutive model

First, it is necessary to specify the assumptions on which the behavior of the material at the grain scale is based. The single crystal constitutive law is assumed to be elastic-plastic, with the plastic deformation only due to crystallographic slip; the other modes of inelastic deformation, such as twinning or phase transformation, are not taken into account in the scope of this study. The local incremental elasto-plastic constitutive law is defined by means of the tangent modulus l relating the nominal stress rate n to the velocity gradient g :

n = l : g (1)
The above velocity gradient is the sum of the total strain rate and spin. The latter two parts can be further divided into an elastic part with the superscript e and a plastic one with superscript p . The plastic parts (plastic strain rate and plastic spin) are related to the slip rates g  by:

where g R and g S denote, respectively, the well-known symmetric and anti-symmetric parts of the Schmid tensor associated with a given slip system g. A rate independent regularization law [START_REF] Franz | Strain localization analysis for single crystals and polycrystals : Towards microstructure-ductility linkage[END_REF] is adopted for the determination of the slip rates:

g g g k    (3) g g g g 0 1 2 gg g c 1 1 1 1 k (1 tanh(k )) (1 tanh(k ( 1))) (1 tanh(k )) H 2 2 2           (4)
where gg H is the self-hardening term and "tanh" stands for the hyperbolic tangent function. Eq. (4) aims at modeling at best a threshold (smooth approximation to the discontinuous step function) without having the numerical problems inherent to the use of sharp functions, and combined with g g g k    provides a regular form of the Schmid law. The expression of the elasto-plastic tangent modulus is then derived as follows (see [START_REF] Lorrain | Critère de ductilité basé sur la perte d"ellipticité du module tangent élastoplastique deduit d'un modèle autocohérent[END_REF][START_REF] Franz | Prédiction de la limite de formabilité des aciers multiphasés par une approche micromécanique[END_REF]

(C R S S )M k R (C )                           (5)
This tangent modulus is composed of elastic and plastic parts, where several convective terms appear due to the large-strain framework. The evolution law of the critical shear stress, in terms of the slip rates, can be expressed as follows (see [START_REF] Lorrain | Critère de ductilité basé sur la perte d"ellipticité du module tangent élastoplastique deduit d'un modèle autocohérent[END_REF][START_REF] Franz | Prédiction de la limite de formabilité des aciers multiphasés par une approche micromécanique[END_REF]):

nslip g gh h c h1 H      (6) with gh gh h c h nslip gk k k1 1 H a ( 2y ) L 2a        ( 7
)
where  is a constant related to the stability of dislocation configurations,  is the shear modulus, and gh a is the anisotropy interaction matrix. To derive the overall response of the polycrystalline aggregate, the self-consistent scaletransition scheme is used (see [START_REF] Lorrain | Critère de ductilité basé sur la perte d"ellipticité du module tangent élastoplastique deduit d'un modèle autocohérent[END_REF][START_REF] Franz | Prédiction de la limite de formabilité des aciers multiphasés par une approche micromécanique[END_REF] for more details). Accordingly, the macroscopic behavior law linking the macroscopic fieldsnominal stress rate N and velocity gradient G -by means of the yet-unknown macroscopic tangent modulus L has the same incremental form as that of the single crystal (Eq. ( 1)):

N = L : G (8)
These macroscopic fields (i.e. velocity gradient and nominal stress rate) are defined as the volume averages of their microscopic counterparts:

vv 11 dv, dv VV   G g N n (9)

Indirect coupling strategy applied to the polycrystalline plasticity model

As discussed before, we adopt a strategy of indirect coupling [START_REF] Kalidindi | An appoximate procedure for predicting the evolution of crystallographic texture in bulk deformation processing of FCC metals[END_REF] between the FE codes Abaqus and LAM3 and the polycrystalline plasticity model, which has been described just above in section 2. This consists in performing the following steps:

 An initial FE calculation that simulates a forming process. In this preliminary simulation, the material parameters are identified with a phenomenological law.  At the end of this preliminary phenomenological FE simulation, the loading history at some appropriately selected material points is recorded.  Then, the above-extracted loading history will serve as prescribed loading for the subsequent simulation, which is based on the polycrystalline model alone (i.e., outside of the FE environment). In what follows, we will denote by "FEC path" the loading path that is extracted from the FE calculation, and by "Linear path" the monotonic or sequential loading path that is applied directly to the polycrystalline model (i.e., without using any preliminary FE calculation).

Applications Indirect coupling: Application to simple loading paths

For the simulation of simple loading paths in the FE code Abaqus, the eight-node reducedintegration C3D8R element is used. Also, the material behavior is taken elasto-plastic with isotropic hardening, with the hardening curve given experimentally. For the polycrystalline model, the identified material parameters are reported in Table 1. where 0  is the initial critical shear stress, 0 g the parameter related to the mean free path of dislocations, c y the critical annihilation distance of dislocations and D the average grain size. In Fig. 1, we compare the polycrystalline model results corresponding to two different ways of prescribing the loading path, namely the loading path extracted from the preliminary FE calculation "FEC path", and the direct application of propositional loading path "Linear path". Then, these polycrystalline results are compared to experimental data and Abaqus calculation.

For the four strain paths reported in Fig. 1, namely uniaxial tension, plane strain tension, simple shear, and Bauschinger (reverse) shear, the results of the polycrystalline model are the same for both definitions of the applied loading path, i.e. "FEC path" and "Linear path". This allowed us to validate our numerical algorithm. The difference between Abaqus results and those of the polycrystalline model is attributable to the material behavior adopted in the FE calculations. In the latter, the hardening is assumed isotropic, while the actual material anisotropy is accounted for in the polycrystalline model along with the associated texture evolution.

Conclusion

In this work, we developed a method of indirect coupling between the FE codes Abaqus and LAM3 and the polycrystalline plasticity model. This method consists in performing a preliminary FE calculation, from which the loading history at some appropriately selected material points is extracted and recorded. The latter will serve to define the loading path to be applied in the subsequent polycrystalline model simulation.

In order to validate this strategy and the associated numerical algorithm, it was first applied to simple loading paths (uniaxial or plane strain tension, simple or Bauschinger (reverse) shear) and then to more complex strain paths that are extracted from skin pass rolling simulations. It is shown that the indirect coupling strategy has several advantages. As compared to a linear loading path method, it can adequately represent the phenomena that occur in a metal sheet under the roll gap.

However, the drawback of this coupling is that it relies on a preliminary FE calculation, in which the material behavior is represented by a phenomenological law. Accordingly, the accuracy of the subsequent micromechanical simulation will depend on the choice of the earlier phenomenological model.

In future works, the texture evolution of the material throughout its loading history will be compared to experimental results.

Fig. 1 .Fig. 3 .Fig. 4 .

 134 Fig. 1. Stress-strain response. (a) Uniaxial tension; (b) Plane strain tension; (c) Simple shear; (d) Bauschinger (reverse) shear.

  ):

	l	ijkl	C	ijkl	ik lj 11 il kj ( ) 22	(	ik lj	il jk	)
			g	g	g	gh h h
			ijpq pq	ip pj	ip pj		mn	mnkl	mn kl

Table 1 :

 1 Identified material parameters for a ferritic steel.

	Parameters	[MPa]		[nm]	[microns] Grains
	Value	100	100	3.25	20	5000

Indirect coupling: Application to rolling skin pass processes

The mesh of the metal sheet and the roll are represented in Fig. 2a, while the strain lines are defined in Fig. 2b. Note that the rolling calculations are performed using the finite element code LAM3 with stationary method. Once again, the material behavior is assumed elasto-plastic with isotropic hardening, with the hardening curve given experimentally. For the polycrystalline model, the identified material parameters are reported in Table 2.

Table 2 :

Identified material parameters for a high strength low alloy steel.

Parameters [MPa] [nm] [microns] Grains

Value 120 90 3.25 20 5000 where 0  is the initial critical shear stress, 0 g the parameter related to the mean free path of dislocations, c y the critical annihilation distance of dislocations and D the average grain size. In Figs. 3 and4, we compare, on a skin pass process, the polycrystalline model results corresponding to the two above-described techniques of loading path application, namely the indirect coupling strategy with the associated loading path extracted from the preliminary FE calculation "FEC path", and the direct application of propositional loading path "Linear path". These polycrystalline results are also compared to the LAM3 FE calculations (see Figs. 34).

From Fig. 3, we can observe that the indirect coupling method along with the associated loading path allows us to follow the strain evolution of the material under roll gap, while the linear loading path method only gives a rough approximation of this evolution. Similar observations can be drawn from Fig. 4, where the indirect coupling method is shown to provide stress evolutions in better agreement with the LAM3 FE results. The difference between the polycrystalline model results and those yielded by the LAM3 FE calculations is also due to the fact that in the latter, the material is assumed isotropic, while in the former, its anisotropy is accounted for along with the associated texture evolution.


A Salahouelhadj
email: abdellah.salahouelhadj@univ-metz.fr@www.univ-metz.fr

H Chalal

F Abed-Meraim

T Balan

ELASTIC-PLASTIC ANALYSES USING THE SOLID-SHELL FINITE ELEMENT SHB8PS AND EVALUATION ON SHEET FORMING APPLICATIONS -COMPLAS XI

Keywords: Solid-Shell Element, Reduced Integration, Physical Stabilization, Assumed Strain Method, Elastic-Plastic Behavior, Sheet Metal Forming

In this contribution, the formulation of the SHB8PS continuum shell finite element is extended to anisotropic elastic-plastic behavior models with combined isotropic-kinematic hardening at large deformations. The resulting element is then implemented into the commercial implicit finite element code Abaqus/Standard via the UEL subroutine. The SHB8PS element is an eight-node, three-dimensional brick with displacements as the only degrees of freedom and a preferential direction called the thickness. A reduced integration scheme is adopted using an arbitrary number of integration points along the thickness direction and only one integration point in the other directions. The hourglass modes due to this reduced integration are controlled using a physical stabilization technique together with an assumed strain method for the elimination of locking. Therefore, the element can be used to model thin structures while providing an accurate description of the various throughthickness phenomena. Its performance is assessed through several applications involving different types of non-linearities: geometric, material and that induced by contact. Particular attention is given to springback prediction for a Numisheet benchmark problem.

INTRODUCTION

During the last decade, considerable effort has been devoted to the development of eightnode solid-shell elements for modeling of thin structures (e.g. [1][2][3][4]). As they use linear interpolation for efficiency reasons, these elements exhibit various locking phenomena which need to be cured in order to preserve the desired accuracy. Nevertheless, compared to conventional shell elements they have many advantages: the use of full three-dimensional constitutive laws, direct calculation of thickness variations, easy treatment to update configurations (no rotational degrees of freedom used), and simple connection with three-dimensional elements since displacements are the only degrees of freedom. For sheet forming applications, key features like double-sided contact and increased accuracy with only one layer of elements through the thickness make these elements particularly attractive.

The reduced integration technique, initiated by the works of Zienkiewicz et al. [5] and Hughes et al. [6], was the first successful solution to alleviate some locking pathologies. Finite elements using this method are very efficient due to their low numerical cost. However, stabilization techniques are needed in order to control the spurious zero-energy deformation modes (or hourglass modes) induced by this reduced integration.

In order to circumvent locking phenomena for three-dimensional low-order elements, several authors have used the enhanced assumed strain (EAS) method, based on Simo and Rifai's pioneer work [7]. The basis of such element formulations is given by the mixed variational principle in which the so-called incompatible strain and stress act as additional independent variables. Recent investigations have combined EAS and reduced integration techniques to derive efficient and accurate elements. As examples, some authors used a fixed number of Gauss points in the thickness direction [1][2][3][4].

The SHB8PS is one such element that has been recently developed [1,2], based on inplane one-point numerical quadrature with eight physical nodes and using an arbitrary number of integration points through the thickness direction. This avoids the use of several layers of elements in order to increase the number of integration points in the thickness, e.g. for metal forming problems. The hourglass modes caused by this reduced integration are efficiently controlled by a physical stabilization technique based on the assumed strain method [8].

In the current contribution, the formulation of the SHB8PS solid-shell finite element is extended to anisotropic elastic-plastic behavior models with combined isotropic-kinematic hardening at large deformations. The resulting element is then implemented into the commercial implicit finite element code Abaqus/Standard via the UEL subroutine. Its good performance is demonstrated through non-linear benchmark problems involving large strains, plasticity and contact. Particular attention is given to springback prediction for a NUMISHEET benchmark problem.

FORMULATON OF THE SHB8PS ELEMENT

Finite element interpolation

SHB8PS is an eight-node, isoparametric hexahedral element with linear interpolation. It has a set of nint integration points chosen along the thickness direction ζ in the local coordinate frame (see Fig. 1).

∑ = = = 8 1 , , , , I I iI I iI i N x N x x ζ η ξ ζ η ξ (1) ( ) ( ) ∑ = = = 8 1 , , , , I I iI I iI i N u N u u ζ η ξ ζ η ξ (2)
Subscript i varies from 1 to 3 and represents the direction of the spatial coordinates. Subscript I varies from 1 to 8.

Discretized gradient operator

First, we introduce the i b (i = 1,..., 3) vectors, representing the derivatives of the shape functions at the origin of the reference coordinate system, defined by Hallquist [START_REF] Hallquist | Theoretical manual for DYNA3D[END_REF] as , (0, 0, 0) 1, 2, 3

T i i b N i = = (3)
The displacement gradient can then be written as follows (see Belytschko and Bindeman [8]):

( ) 4 , , , 1
T T T T i j j j i j j i u b h d b h d α α α α α γ γ = = + ⋅ = + ⋅       ∑ (4)
where i d are the nodal displacement vectors. The functions h α and vectors α γ (α =1,...,4) are given by

1 2 3 4 , , , h h h h ηζ ζξ ξη ξηζ = = = = (5) 1 1 ( ) 8 T j j j h h x b α α α γ = = - ⋅       ∑ (6)
The discretized gradient operator can be written as

0 0 0 0 0 0 T T x x T T y y T T z z T T T T y y x x T T T T z z y y T T T T z z x x b h b h b h B b h b h b h b h b h b h α α α α α α α α α α α α α α α α α α γ γ γ γ γ γ γ γ γ + + + = + + + + + +                     , , , , , , , , , 0 0 0

Stabilization and assumed strain method

The particular location of the integration points along a line generates six so-called hourglass modes. The control of the hourglass modes of the SHB8PS element is achieved by adding a stabilization component STAB K to the element stiffness matrix e K . This part is drawn from the work of Belytschko and Bindeman [8], who applied an efficient stabilization technique together with an assumed strain method. The stabilization forces are consistently derived in the same way. Moreover, the discretized gradient operator is projected onto an appropriate sub-space in order to eliminate shear and membrane locking. In this approach, the i b vectors (Eq. ( 3)) are replaced by the mean value of the derivatives of the shape functions over the element, denoted by ˆi b , as proposed by Flanagan and Belytschko [START_REF] Flanagan | [END_REF]: B part defined identically but where α varies only from 1 to 2. Then, the remaining part 34 B of B ˆ, which vanishes at the integration points, is further projected as 34 B .

One can project the B ˆ operator onto a B ˆ operator as:

12 34 ˆˆB B B = + (9)
where 34 B is given by:

0 0 ˆ0 0 0 0 0 0 0 0 T x T y T z T x h h h B h α α α α α α γ γ γ γ = = =                   ∑ ∑ (10)
The stiffness matrix e K takes the form:

12

ˆˆ e T e Geom STAB Geom ep K B C B d K K K K Ω = ⋅ ⋅ Ω + = + + ∫ (11)
where the first term 12 K is evaluated at the integration points as

Ω Ω Ω = ⋅ ⋅ Ω + ⋅ ⋅ Ω + ⋅ ⋅ Ω ∫ ∫ ∫ (13)
In a similar way, the internal forces of the element can be written as

where STAB f represents the stabilization forces. The stabilization terms are calculated in a co-rotational coordinate system [8].

NUMERICAL EXAMPLES

Pinched cylinder

In this example, the elastic-plastic deformation of a cylinder subjected to two opposite concentrated loads in the middle of the structure and bounded by rigid diaphragms on its extremities is considered. This problem has been investigated by a number of authors like [2,[START_REF] Simo | [END_REF][13][START_REF] Sori\'{c} | {a}tzig[END_REF][START_REF] Brank | Damjani\'{c}, FB[END_REF][16].

The undeformed mesh and boundary conditions are shown in Figure 2. The geometry is characterized by the length L=600 mm, the radius R =300 mm and thickness t = 3 mm. Due to symmetry, only one eighth of the cylinder is discretized. At the ends of the cylinder, the rigid diaphragms prevent any displacement in the radial directions. Material properties are the elasticity modulus E = 3000 MPa, Poisson's coefficient ν = 0.3 and initial yield stress σ 0 = 24.3 MPa. A linear isotropic hardening law is adopted and can be written as:

P eq Y Hε σ σ + = 0 (15)
Where P eq ε is the equivalent plastic strain and H the linear hardening coefficient taken equal to 300 MPa. Figure 3 shows the load versus the vertical deflection. For the same mesh (40×40×1 elements), the results of the SHB8PS solid-shell element are compared to those of the SC8R solid-shell element and the S4R shell element, along with the reference solutions obtained by [13] and [17]. The calculations using the S4R element failed at a certain loading level, while the SC8R element is too stiff in this test problem and also converges more slowly. For this particularly discriminating test, the curve using SHB8PS element is clearly in better agreement with the 3D investigations of [13] and [17].

Unconstrained cylindrical bending

The example of the unconstrained cylindrical bending test proposed as springback benchmark in NUMISHEET 2002 is studied [START_REF] Meinders | 6th international conference and workshop on numerical simulation of 3D sheet metal forming processes[END_REF]. This application allows us to evaluate the performance of the SHB8PS element, implemented in Abaqus/Standard, in presence of geometric, material and contact non-linearities. This benchmark involves a bendingdominated deformation since there is no blank holder. The problem has complex contact boundary conditions during the forming process and the springback after forming is severe. The geometry of the problem is illustrated in Figure 4 and the geometric parameters are summarized in Table 1.

The material under investigation is a High Strength Steel, which is supposed elastic-plastic with isotropic hardening following Swift law: The amount of springback is quantified by the angle θ as defined in Figure 5. This angle is measured after forming at the maximum punch displacement and after springback. The tools are defined as analytical rigid surfaces. The SHB8PS element is compared with both solid and shell elements. Indeed, it is wellknown that in applications of sheet metal forming, shell elements have difficulties in dealing with double-sided contact -while conventional solid elements require several element layers to capture bending effects. In the present work, the simulations carried out with the SHB8PS element use only one element layer through the thickness. For symmetry reasons, only one quarter of the blank is discretized by means of 150 SHB8PS elements in the length and only one element over the width of the sheet. The analysis with the SHB8PS element is carried out using five Gauss points in the thickness direction because elastic-plastic applications require, in general, five integration points in minimum to describe the strongly non-linear throughthickness stress distribution [2].

( ) 0 Y n P eq K σ ε ε = + (16)
In order to validate the proposed solid-shell element, its predictions are compared to the experimental results of the Numisheet 2002 benchmarks. Two elements from the element library of the Abaqus code are also used in the comparison: the shell element S4R and the 3D continuum element C3D8I. Again, 150 uniformly distributed elements are used in the length direction for these two elements. However, two layers of C3D8I elements are required in the thickness direction in order to represent the stress distribution due to bending with sufficient accuracy. Also, ten C3D8I elements are used along the width direction in order to keep their aspect ratio in acceptable limits. Figure 4 displays the punch force versus punch displacement curves predicted by the three elements, along with the experimental results (BE-1 to BE-4) from Meinders et al. [START_REF] Meinders | 6th international conference and workshop on numerical simulation of 3D sheet metal forming processes[END_REF].

Figure 6 shows that the numerical results obtained with SHB8PS element are the closest to the experimental results and they lay close to the solid element predictions. The slight differences between the two may be due to the different number and distribution of integration points along the thickness direction. The S4R element has too soft behavior with respect to SHB8PS and C3D8I elements.

The springback angles are also investigated, as they were also experimentally measured [START_REF] Meinders | 6th international conference and workshop on numerical simulation of 3D sheet metal forming processes[END_REF]. The springback phenomenon is particularly exacerbated in this unconstrained bending application, as illustrated in Figure 7. Table 2 summarizes the opening angles before and after springback for elements SHB8PS, C3D8I and S4R, compared to experiments. The simulated values with SHB8PS and C3D8I elements are close to each other and the closest to experiments. Comparing the numerical results to the experimental ones, the good performance of the SHB8PS solid-shell element is confirmed.

CONCLUSIONS

An extended version of the solid--shell finite element SHB8PS has been implemented into the implicit finite element code Abaqus/Standard via the UEL subroutine. The formulation of this element employs a combination of the reduced integration scheme with the assumed strain method and a specific projection to eliminate locking phenomena. The resulting hourglass modes are controlled using a physical self-adapting stabilization procedure. This version of the SHB8PS element can deal with problems involving anisotropic elastic-plastic behavior at large deformations and double-sided contact between sheet and tools, which are typical in sheet metal forming applications.

Figure 1 :

 1 Figure 1: SHP8PS reference geometry

  α γ are replaced by vectors ˆα γ where the i b vectors are simply substituted by ˆi b . A modified discretized gradient operator B ˆ can be constructed in the same way. It can be shown that the terms of the B ˆ operator vanish for α = 3,4. In other words, the B ˆ operator reduces to its 12

  the Jacobian of the transformation between the reference and the current configurations; ( ) I ω ζ is the corresponding weight, while ep C modulus. The geometric stiffness matrix Geom K is due to the non-linear (quadratic) part of the strain tensor and STAB K represents the stabilization stiffness given by equation:

Figure 2 :

 2 Figure 2: Geometry and boundary conditions for the pinched cylinder

Figure 3 :

 3 Figure 3: Load deflection curves for elastic-plastic pinched cylinder

  the hardening parameters of the material. The Young modulus E = 2.175 × 10 5 MPa and the Poisson ration ν = 0.3. Further K = 645.24 MPa, n = 0.25177 and 0 ε = 0.0102. The friction coefficient of the interaction between surfaces punch-sheet and diesheet is µ = 0.14812.

Figure 4 : 1 :

 41 Figure 4: Tool geometry for the unconstrained bending problem

Figure 5 :

 5 Figure 5: Definition of the angle to measure springback for the unconstrained cylindrical bending problem

Figure 6 :

 6 Figure 6: Punch force vs. punch displacement curves for High Strength Steel

Figure 7 :

 7 Figure 7: Deformed shape of the sheet in the unconstrained bending problem

Table 2 :

 2 Measured and simulated opening angles before and after springback

			Experimental			Simulated
		BE-01	BE-02 BE-03	BE-04	SHB8PS C3D8I	S4R
	Forming	22.7707	22.0064 23.0255 20.8599 23.0692	22.5820 33.3078
	Springback 37.4212	35.6787 30.9036 35.3636 36.3952	32.0832 43.9071

The performance of SHB8PS element has been shown through two numerical examples involving various types of non-linearities: geometric, material and contact. Indeed, at equivalent mesh density, SHB8PS performs at least as well as the most accurate (and expensive) solid elements. However, this accuracy is achieved at a lower cost by using one single layer of SHB8PS elements and simply adjusting the number of integration points through the thickness. This feature makes the element very competitive for sheet metal forming.