Measurement of the Y(1S) pair production cross section and search for resonances decaying to Y(1S)mu(+)mu(-) in proton-proton collisions at root s=13 TeV

(1)

Contents lists available atScienceDirect

Physics

Letters

B

www.elsevier.com/locate/physletb

Measurement

of

the

Y(1S)

pair

production

cross

section

and

search

for

resonances

decaying

to

Y(1S)μ

+

μ

−

in

proton-proton

collisions

at

√

s

=

13 TeV

.

The

CMS

Collaboration

CERN,Switzerland

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory: Received15February2020 Receivedinrevisedform2June2020 Accepted22June2020

Availableonline26June2020 Editor: M.Doser

Keywords: CMS Upsilon

ThefiducialcrosssectionforY(1S)pairproductioninproton-protoncollisionsatacenter-of-massenergy of13 TeV intheregionwherebothY(1S)mesonshaveanabsoluterapiditybelow2.0ismeasuredtobe 79±11(stat)±6(syst)±3(B)pb assumingthemesonsareproducedunpolarized.Thelastuncertainty corresponds to the uncertainty inthe Y(1S) meson dimuon branching fraction. The measurement is performedinthe finalstatewith fourmuonsusingproton-protoncollisiondata collectedin2016by the CMS experimentatthe LHC, corresponding toan integratedluminosity of35.9 fb−1.Thisprocess servesasastandardmodelreferenceinasearchfornarrowresonancesdecayingtoY(1S)μ+μ− inthe samefinalstate.Sucharesonancecouldindicatetheexistenceofatetraquarkthatisaboundstateof twob quarksandtwob antiquarks.¯ Thetetraquarksearchisperformedformassesinthevicinityoffour timesthebottomquarkmass,between17.5and19.5 GeV,whileagenericsearchforotherresonancesis performedformassesbetween16.5and27 GeV.Nosignificantexcessofeventscompatiblewithanarrow resonanceisobservedinthedata.Limitsontheproductioncrosssectiontimesbranchingfractiontofour muonsviaanintermediateY(1S)resonancearesetasafunctionoftheresonancemass.

1. Introduction

Quarkonium pair production is an important probe of both perturbative and nonperturbative processes in quantum chromo-dynamics.Experimental studies ofthisprocess can provide valu-able information about the underlying mechanisms of particle production andimprove our understanding of numerous physics processes that are treated asbackgrounds in searches and mea-surements. Quarkonium pairs may originate from single-parton scattering(SPS) ordouble-parton scattering(DPS). These produc-tion mechanisms can be separated experimentally since the DPS production is characterized, among other features, by more for-ward and separated mesons. The analysis of nonperturbative ef-fects is easier for quarkonium states composed of b quarks, as their large masses allow them to be approximated as nonrela-tivisticsystems [1]. The CMS Collaboration observed forthe ﬁrst timetheproductionofapairofY

(

1S

)

mesons,usingproton-proton datacollected ata center-of-massenergyof8 TeV [2].ThisLetter presentsa measurement of the Y

(

1S

)

pair production cross sec-tion at a center-of-mass energy of 13 TeV. The cross section is measuredintheﬁducialregionwherebothY

(

1S

)

mesonshavean

_E-mail_address:_cms_-publication_-committee_-chair_@cern_.ch.

absoluterapiditybelow2.0,usingtheﬁnalstatewithfourmuons. Additionally, theDPS contributiontothe process ismeasured for theﬁrsttime.

TheY

(

1S

)

pairproductioncan serveasareferenceinsearches fortetraquarks orgeneric resonances withmassesclose totwice theY

(

1S

)

mesonmass.AlightresonancedecayingtoaY

(

1S

)

me-son anda pairof leptons mightbe the signature ofa tetraquark characterized as a bound state of two b quarks and two b anti-

¯

quarks,especially ifits mass isbelow twice the

η

b mass [3–13].

In this Letter, in addition to the measurement of the Y

(

1S

)

pair productioncrosssection,wedescribeasearchfortetraquarkswith massesbetween17.5and19.5 GeV,sincebbb

¯

b tetraquarks

¯

would beexpectedtohaveamassaroundfourtimesthatofthebottom quark.Agenericsearchfornarrowresonanceswithmassbetween 16.5 and 27 GeV and decaying to a Y

(

1S

)

meson and a pair of muonsisalsopresented.Theﬁnalstateisthesameasforthe mea-surementoftheY

(

1S

)

pairproductioncrosssection,andasimilar eventselectionisused.TheY

(

1S

)

pairproductionisabackground totheresonancesearch.

The LHCb Collaboration searched for bbb

¯

b tetraquarks

¯

using datacollectedatcenter-of-massenergiesof7,8,and13 TeV, with-out ﬁndinganyhintofasignal [14]. Thisanalysisprobesa

kine-https://doi.org/10.1016/j.physletb.2020.135578

(2)

maticregionthatisnotaccessiblewiththeLHCbdetectorand ex-tendsthecoveredmassrangeinthecontextofthegenericsearch. The Y

(

1S

)

pair production ﬁducialcross section measurement andtheresonancesearcharebasedonproton-protoncollisiondata collectedin2016atacenter-of-massenergyof13 TeV bytheCMS experimentattheCERNLHC,correspondingtoanintegrated lumi-nosityof35.9 fb−1.

2. TheCMSdetector

The central feature of the CMS apparatus is a superconduct-ing solenoidof 6 m internal diameter, providinga magnetic field of3.8 T.Withinthesolenoidvolume,thereareasiliconpixeland strip tracker,a leadtungstatecrystalelectromagnetic calorimeter, andabrass andscintillatorhadron calorimeter,eachcomposedof abarrelandtwoendcapsections.Forwardcalorimetersextendthe pseudorapiditycoverageprovidedbythebarrelandendcap detec-tors.Muonsaredetectedingas-ionizationchambersembeddedin thesteelflux-return yokeoutsidethesolenoid. Events ofinterest are selected using a two-tiered trigger system [15]. A more de-tailed descriptionof theCMSdetector, together witha definition ofthecoordinatesystemusedandtherelevantkinematicvariables, canbefoundinRef. [16].

Muons are measured in the range

|

η

| <

2

.

4, with detection planes made using three technologies: drift tubes, cathode strip chambers,andresistiveplatechambers.Matchingmuonstotracks measured inthesilicontrackerresultsin arelative pT resolution

intherange0.8–3.0%formuonswith

p

T lessthan10 GeV [17]. 3. Simulatedsamples

TheY

(

1S

)

pairproductionsignalissimulatedusingthe pythia 8.226generator [18],separately fortheSPS andDPSmechanisms, undertheassumptionthat themesonsareproduced unpolarized. The DPSsample isproduced by generatingtwo hard interactions withcolor-singletproductionofbottomoniumstates viagg

→

bb

¯

orcolor-octetproductionofbottomoniumstatesviaqq

→

bb.

¯

The invariant massdistributionof themesonpairandofthe rapidity separationbetweenthemesonsareusedtoextractthefractionof DPSproduction,asdetailedinSection5.Forthismeasurement,the distributionsofthesevariablesfortheSPSprocessaretakenfrom the next-to-leading-order (NLO*) calculation with a cutoff color-singletmechanism(CSM) [19–21] using HELAC-Onia 2.0.1 [22,23].

The signal of a narrow resonance decaying to a Y

(

1S

)

meson anda pair ofmuonsis modeled using differentphysics assump-tionsdependingonthenatureoftheresonance:

•

a bottomonium state withthe properties of the

χ

b1(1P),

as-suminga phase-space decayto a Y

(

1S

)

mesonanda pair of muons,usingthe pythia 8.226generator;

•

a scalarparticleproduced in gluonfusion, usingthe JHUGen generator [24–27];

•

a pseudoscalar particle produced in gluon fusion, using the JHUGengenerator;

•

aspin-2particle producedingluon fusion,usingthe JHUGen generator.

Thesignalsaregeneratedassumingthenarrow-width approxima-tion. The

χ

b1(1P)sample is used to modelthe tetraquark signal,

forwhichnodedicatedgeneratorexists.Theothersamples corre-spondtothesignalsinthegenericsearch overanextendedmass range.Foreachmodel,fourresonancemassvaluesaresimulated: 14,18, 22, and26 GeV. Since the signal acceptance fallssteeply around andbelow 14 GeV in the simulated samples,the probed mass range in this analysis is restricted to stay well above this

mass threshold.The differentmasspoints areused tointerpolate andextrapolatethesignalmodeloverthewholemassrange.

The pythia generatorwiththetuneCUETP8M1 [28] is usedto modelthepartonshower andhadronizationprocesses.Generated events are processed through a simulation of the CMS detector basedon Geant4 [29].

4. Eventselectioncriteria

The event reconstruction is based on the particle-ﬂow algo-rithm [30], which identiﬁes individual particle candidates using informationfromalltheindividualsubdetectors.Muonsare recon-structedbycombininginformationfromthesilicontrackerandthe muonsystem [17].

Events are selected with a trigger that requires the presence ofthreemuons.Among thesemuons,two musthaveaninvariant mass compatible with a Y resonance (8

.

5

<

m2μ

<

11

.

4 GeV) at

triggerlevel,andthedimuonvertexﬁtprobability,calculatedusing the

χ

2 _and_the_number_of_degrees_of_freedom_of_the_ﬁt,_must_be

greaterthan0.5%.

Oﬄine, we require each event to have four reconstructed muonswith pT

>

2 GeV and

|

η

| <

2

.

4.Thesemuonsarerequired

to satisfy the globalor particle-flow muon identification criteria described inRef. [17]. About 25% of simulated signal eventsand about30% of dataeventshave morethan foursuch muons. Pos-siblecombinationsoffourmuontracksare refitwithaconstraint tocomefromacommonvertex,andthe

χ

2 _probability_of_the_ﬁt

isdetermined.Thecombinationoffourmuonswiththelargest

χ

2

probability ischosen.Forsimulatedsignaleventswithmorethan fourreconstructedmuons, thecorrectmuonsarechosen inabout 98%ofcases.Amongthefourmuons,atleastthreeneedtobe as-sociated withthe trigger-level objects. At leasttwo muons must beassociatedwiththeobjectsthatpassedtheY mass compatibil-ityandvertexcriteriaofthetrigger,andtheyarepairedtogether. If thereare more thantwo such muons, which happensfor2 to 35% ofsimulatedsignaleventsdependingontheresonancemass, those thathaveopposite-sign(OS)chargesandaninvariant mass closesttotheworld-averageY

(

1S

)

mass [31] arepairedtogether.

After selectingthe bestcombinationoffourmuonswith pT

>

2 GeV,the

p

Tthresholdisraisedto2.5 GeV fortheselectedmuons.

Theﬁnalselectionrequiring

p

T

>

2

.

5 GeV reducesthebackground

from misidentified muons by about a factor of two. The muons arerequiredtosatisfythemediummuonidentificationcriteria de-scribed in Ref. [17]. Both pairs of muons have to be composed of OSmuons. The vertexfit

χ

2 _probability _of_the _four _muons_is

required to be greater than 5%, whereas that of the Y

(

1S

)

can-didate is required to be above 0.5%, similar to the requirement already imposed at trigger level. The muons are required to be separated from each other by at least

R

=

√

(

η

)

2

+ (φ)

2

=

0

.

02,where

η

and

φ

arethedifferencesinpseudorapidityand azimuthal angle between the muons. The positively (negatively) chargedmuonfromoneofthepairscanbepairedwiththe nega-tively(positively)chargedmuonoftheotherpairtoformso-called alternative pairs of OS muons. If one of these alternative pairs has an invariantmass compatiblewitha J

/

ψ particle within two standard deviations of theexperimental resolution, which ranges between about 0.03 and 0.12 GeV depending on the muon pair kinematics, the event is discarded from the analysis. Events are alsodiscardediftheycontaintwo OSpairs ofmuonswith invari-antmasslessthan4 GeV.

Theselection criteriadetailedabovearecommonforthe mea-surementoftheY

(

1S

)

pairproductioncrosssectionandthesearch for a resonant signal. The criteria that differ between the mea-surement and the search are described in the following. In the measurement ofthe Y

(

1S

)

pair ﬁducial cross section, the recon-structed absolute rapidity of both muon pairs is required to be

(3)

lessthan2.0.Inaddition,formuonswith

|

η

| <

0

.

9,the

p

T

thresh-old israised to 3.5 GeV. Centralmuons withtransverse momen-tumbelow3.5 GeV havea highprobability of beingabsorbedin the calorimeter orundergoing signiﬁcant multiple scattering be-forereaching themuon detectors.Thisselection criterionreduces thesystematicuncertainty inthe muon reconstruction relatedto the detectorsimulation. It is, however, not used in the resonant searchbecauseitwouldstronglyreducethesignal acceptancefor thelower-mass signal range.In the resonancesearch, the invari-ant mass of the Y

(

1S

)

candidate is required to be within two standarddeviationsoftheexperimental resolutionfromtheY

(

1S

)

mass [31], where the resolution varies between about 0.06 and 0.15 GeV dependingontheevent.

The mass range of interest is known a priori for the search of a bbb

¯

b tetraquark

¯

signal. In this case, all the selection crite-riadescribedabove havebeendeterminedandﬁxedin ablinded way,usingsimulationandwithoutlookingatdataeventswithfour muonshavinganinvariantmassbetween17.5and19.5 GeV.

5. MeasurementoftheY

(

1S

)

pairproductioncrosssection

The methodology used to measure the Y

(

1S

)

pair production crosssectionisdetailedinSection5.1.Afterdiscussingthe system-aticuncertaintiesinSection5.2,theresultsofthemeasurementof theinclusiveY

(

1S

)

pairproductionﬁducialcross sectionare pre-sented in Section 5.3. Nonisotropic decays of the Y

(

1S

)

mesons would change the measured cross section. Section 5.4 describes how the crosssection wouldvary for nonzero values ofthe po-larization parameters. Finally, the DPS and SPS mechanisms can beseparatedexperimentallybymeasuring theY

(

1S

)

pair produc-tion cross section in bins of the rapidity difference betweenthe mesons,

|

y

(

Y

(

1S

),

Y

(

1S

))

|

, andofthe invariantmass ofthe me-sonpairs,

m

Y(1S)Y(1S).AmeasurementoftheDPS-to-inclusivecross

sectionratiointheﬁducialregionispresentedinSection5.5.

5.1. Methodology

TheY

(

1S

)

pairproductioncrosssectionismeasuredinthe fidu-cialregionwherebothmesonshaveanabsoluterapiditybelow2.0. Noother requirementisappliedto definethefiducialregion.The fiducialcrosssection,

σ

ﬁd,canbeexpressedas:

σ

ﬁd

=

Ncorr

LB

2

,

(1)

where

N

corr_is_the_number_of_signal_events_corrected_for_the

accep-tanceandeﬃciencyoftheselection,

L

istheintegratedluminosity, and

B

stands for

B(

Y

(

1S

)

→

μ+μ−

)

= (

2

.

48

±

0

.

05

)

% [31],which is the branching fraction of the Y

(

1S

)

meson decay to a pair of muons. To extract Ncorr _from_the _data, _we _perform _an _extended

unbinnedtwo-dimensional(2D)maximumlikelihoodﬁtofthe in-variantmassdistributionsoftwoOSmuonpairs,whereallevents are weighted for the acceptance and eﬃciency on an event-by-eventbasisbytheweight

ω

,deﬁnedas:

ω

=

A1A2

1reco

2reco

1

− (

1

−

₁vtx

)(

1

−

₂vtx

)

evt

−1

,

(2)

wherethedifferenttermsaredescribedbelow:

•

A, theprobability fora Y

(

1S

)

mesonwithan absolute rapid-ity below 2.0 anddecaying to a pair of muons to have two muons in the geometrical acceptanceof the detector(muon

|

η

| <

2

.

4);No strongcorrelationbetweentheacceptance val-uesofthetwomesonsarefoundwithaclosuretestdescribed inSection

5.2

,andthetotalacceptanceisthereforecomputed astheproductoftheper-mesonweights;

•

reco_,_the _probability _for_a _Y

₍

_1S

₎

_meson_with_an _absolute

ra-pidity below2.0 anddecayingto a pairof muonseach with

|

η

| <

2

.

4 tohavetworeconstructedmuonspassingthe identi-ﬁcationandkinematiccriterialistedinSection4;

•

vtx_,_the_probability_for_a_Y

₍

_1S

₎

_meson_passing_the_acceptance

reconstructioncriteriaoutlinedinitems2and3tohavea ver-texﬁt

χ

2 _probability_above_0.5%;

•

evt_,_the_probability_for_an_event_where_both_Y

₍

_1S

₎

_candidates

passallthecriteriaofitems2and3,andatleastoneofthem passesthevertexﬁt

χ

2_probability_criterion_of_item_4,_to_pass

thefollowingevent-levelcriteria:thetriggerrequirements,the four-muonvertexﬁt

χ

2_probability_above_5%,_and_the_absence

ofOS dimuonpairs withan invariant masswithin two stan-darddeviationsoftheworld-averageJ

/

ψ mesonmass [31]. Theﬁrstthreeitemsintheabovelistarecalculatedasafunction oftheY

(

1S

)

rapidityand

p

T.Thevaluesof A, reco,and

vtx,range

between0.47and1.00, 0.23and0.88, and0.81and0.98, respec-tively,depending onthe Y

(

1S

)

rapidity and pT.The factor

evt is

calculated asa function of the pT of both Y

(

1S

)

candidates, and

rangesbetween0.33and0.65. ThesubscriptindicesinEq. (2) in-dicate theY

(

1S

)

candidate towhich theweight corresponds.The factor

vtx_enters_the_formula_differently_from_the_other_acceptance

andeﬃciencytermsbecausethedimuonvertexﬁt

χ

2 _probability

criterion needs to be satisﬁed by at least one of the two Y

(

1S

)

candidates, but not necessarily by both. The weight

ω

is com-puted onan event-by-eventbasis, usingthe kinematicquantities ofthereconstructedY

(

1S

)

candidatesindata.Theyareestimated fromsimulationasefficiencymapsandaresimilarfortheSPSand DPSproductionmodes,despitedifferentcorrelationsbetweenthe mesons. Data-to-simulation correctionsfor the trigger andmuon identification efficiencies are takeninto account inthe computa-tionofNcorr_.

Inabout3%ofcases,thefourreconstructedmuonsarenot cor-rectly paired in the SPS and DPS Y

(

1S

)

pair simulations. These events cannot be identified as part of the signal by the 2D fit since theirdistribution issimilarto thatofthe floating combina-torialbackground.Therefore,thevalue

N

corr_extracted_from_the_ﬁt

iscorrectedby

+

3% totakeintoaccountthesemispairings. Inthe 2D ﬁt,themuonsare paired asdescribed inSection 4, andthe invariant massesofthe two pairs are randomlydenoted

m12 and m34. The signal model corresponds to Y

(

1S

)

+

Y

(

1S

)

events,whereasthebackgroundmodelisthesumofthefollowing physicsprocesses:

•

Y

(

2S

)

+

Y

(

2S

)

;

•

Y

(

3S

)

+

Y

(

3S

)

;

•

Y

(

2S

)

+

Y

(

1S

)

;

•

Y

(

3S

)

+

Y

(

1S

)

;

•

Y

(

1S

)

+

combinatorial background;

•

Y

(

2S

)

+

•

Y

(

3S

)

+

•

combinatorialbackground

+

combinatorialbackground. TheshapeoftheinvariantmassdistributionfortheY

(

1S

)

com-ponent isdetermined froma 2D ﬁt ofthe two dimuoninvariant massesintheY

(

1S

)

pairSPSsimulation.Theresultsareveriﬁedto be compatiblewiththose ofa ﬁt performedusingthe simulated DPSevents,evenifthemuonrapiditydistributionsdifferbetween productionmodes. The m12 and

m

34 distributions are ﬁttedwith

the sum of two same-mean Crystal Ball functions, which corre-spondto a power lawtailadded to a Gaussian core. Thisallows the radiative tails ofthe distributions to be well modeled. Fig. 1

showstheprojection ofthe2D ﬁtonthe

m

12axisforY

(

1S

)

Y

(

1S

)

simulated events. The projection on the m34 axis is statistically

(4)

Crys-Fig. 1. Projectionofthe2Dﬁt(line)tothem12invariantmassdistribution(points)

fortheSPSY(1S)Y(1S)simulation.Theverticalbarsonthepointsshowthe statisti-caluncertaintyonly.ThemassdistributionismodeledwiththesumoftwoCrystal Ballfunctionswiththesamemean.

talBallfunctionsinsimulationiscompatiblewithin onestandard deviationwiththeworld-averagemassoftheY

(

1S

)

meson,while thefull widthathalf maximumisabout0.19 GeV, whichis sev-eralorders of magnitudelarger than the world-average width of theY

(

1S

)

meson [31] becauseofthelimiteddetectorresolution.

ThecontributionsfromY

(

2S

)

andY

(

3S

)

mesonsare small,and the dimuon invariant mass distributions for these mesons are takenfromacontrol regionindatawitheventswithtwo muons andtwo additionaltracksthatdonot correspondtomuon candi-dates.BothprocessesaremodeledwithaGaussianfunction.

Thecombinatorialbackgroundcomponentsinthe

m

12and

m

34

distributionsaremodeledwithsecond-orderChebychev polynomi-alswithidenticalparameters. The numberofdegreesoffreedom hasbeen determined with a Fisher F-test [32], where the distri-butionofthe combinatorialbackgroundisfoundby invertingthe muonpairassociationinthesignal region.Theparametersofthe polynomialare freetoﬂoatinthe2D ﬁttodatainthesignal re-gion,detailedinSection5.3.

In the 2D ﬁt to the data performed in the signal region, the freeparametersarethenormalizationsofalltheprocessesandthe parametersofthecombinatorialbackgroundmassdistribution.The function parameters ofthe Y

(

1S

)

, Y

(

2S

)

,andY

(

3S

)

signal shapes areconstrainedwithintheiruncertainties.

5.2. Systematic uncertainties

The normalization uncertainties that affect the measurement arethefollowing:

•

2.5%uncertaintyintheintegratedluminosityforthe2016 run-ningperiod [33],whichappearsinEq. (1).

•

0.5% uncertainty per muon in the eﬃciency of the muon identiﬁcation and tracking, measured with a tag-and-probe method [17]. It sumsup to2% per eventbecause the uncer-taintiesareassumedtobecorrelatedforthefourmuonssince theymostly originate fromthesamesource.Thisuncertainty isrelatedtotheterm

reco_in_the_weight

_ω

_.

•

1% uncertainty in the vertexﬁt

χ

2 _probability _criterion,

de-termined by comparing background-subtracted observed and simulated distributions of the vertex ﬁt

χ

2 _probability _for

events witha Y

(

1S

)

meson andtwo nearby tracks. This un-certaintyisrelatedtotheterm

vtx_in_the_weight

_ω

_.

•

2% uncertainty per muon matched to trigger objects in the trigger eﬃciency, measured with a tag-and probe method,

Table 1

SystematicuncertaintiesconsideredintheY(1S)pairproductioncrosssection mea-surement.Thelastcolumngivestheassociatedabsoluteuncertaintyinthe mea-surementof

σ

ﬁd.

Uncertainty source Uncertainty (%) Impact onσﬁd( pb)

Integrated luminosity 2.5 2.0

Muon identiﬁcation 2.0 1.6

Trigger 6.0 4.7

Vertex probability 1.0 0.8

B(Y(1S)→μ+μ−) 4.0 3.2 Signal and background models 1.2 1.0

Method closure 1.5 1.2

Total 8.1 6.4

summing up to 6% per event because the uncertainties are assumedtobecorrelatedforthethreemuonsrequiredat trig-ger level. Thisuncertainty is related to the term

evt _in _the

weight

ω

.

Thesenormalizationuncertaintiespropagatedirectlyintoidentical uncertainties intheY

(

1S

)

pairproductioncrosssection. Addition-ally,theuncertaintyof2%inthe

B(

Y

(

1S

)

→

μ+μ−

)

branching frac-tion,whichisusedtocompute

N

corr _based_on_{Eq. (}₁_),_results_in_a

4%uncertaintyintheY

(

1S

)

pairproductioncrosssection measure-ment.

The parameters of the combinatorial background are freely ﬂoating, whiletheparameters oftheY

(

1S

)

Y

(

1S

)

distributions are constrainedwithintheuncertaintiesobtainedfromtheﬁtto simu-latedevents.Anuncertaintyof0.2%inthemuonmomentumscale is propagatedasan uncertaintyin themeanof theY

(

1S

)

model. These uncertainties inthe signal andbackgroundmodeltogether contributeanuncertaintyof1.5%intheY

(

1S

)

pairproductioncross sectionmeasurement.

Theconsistencyofthemethodtoobtain

N

corr_is_checked_by

ap-plyingtheeﬃciencyandacceptanceweightstotheeventsselected in simulation,andcomparing thecomputed Ncorr _to _the_number

ofeventsgeneratedintheﬁducialregion beforeapplyingany se-lectioncriterion. Thistestisperformedforboth theSPSandDPS simulations usingthe correction mapsderived fromone sample, theotherone,ortheircombination.Usingthecombinedmap,the weighted DPSyieldhasadeviationof

(

−

1

.

3

±

3

.

7

)

% withrespect tothegeneratedyield,andthecorrespondingdeviationfortheSPS sample is

(

−

0

.

6

±

1

.

5

)

%. Thelevelofclosureissimilarlygoodfor bothproductionmodesdespiteaverageeventweightsdifferingby morethana factorof3becauseofthekinematicdifferences. The weighted numberof dataeventsused tocompute the Y

(

1S

)

pair production crosssection isincreasedby 1% toallow fora poten-tialnonclosure,andanuncertaintyof1.5%isassociated withthis correction.

ThesystematicuncertaintiesaresummarizedinTable1.

5.3. Measurement of the ﬁducial cross section

The 2D unbinned ﬁt to the m12 vs. m34 distribution yields Ncorr

=

1740

±

240 fortheY

(

1S

)

Y

(

1S

)

process.Theprojectionson both dimensionswithall theﬁtcomponentsare showninFig.2. This number of events can be translated into an inclusive cross section forthe Y

(

1S

)

Y

(

1S

)

process in the ﬁducial region deﬁned such thatbothY

(

1S

)

mesonshaveanabsoluterapiditybelow2.0. Takingintoaccountthestatisticalandsystematicuncertainties de-scribed in Section 5.2, and assuming unpolarized Y

(

1S

)

mesons, theinclusiveﬁducialcrosssectionismeasuredtobe:

σ

ﬁd

=

79

±

11 (stat)

±

6 (syst)

±

3

(

B

)

pb

,

(3)

where the last uncertainty comes from the uncertainty in the Y

(

1S

)

dimuonbranchingfraction.

(5)

Fig. 2. Thetwoprojectionsandthe resultofthe2Dfittothemuonpair invari-antmasses.Eacheventiscorrectedforacceptanceandefficiency.TheY(1S)pair productionsignalisshownasafilledarea.Thecontributionsfromthe combinato-rialbackground,andfromeventswithaY(1S)mesonandapairofcombinatorial muons,withaY(2S)mesonandtworeconstructedmuonsfromanyorigin,and withaY(3S)mesonandtworeconstructedmuonsfromanyorigin,areoverlaid.

TheCMSCollaborationpreviouslymeasured,inthesame ﬁdu-cialregion,theY

(

1S

)

Y

(

1S

)

productioncrosssectionata center-of-massenergyof8 TeV to be 69

±

13(stat)

±

7(syst)

±

3

(B)

pb [2]. Assumingall uncertainties areuncorrelatedwiththose inthe re-sultpresented inthisLetter exceptthat inthebranchingfraction oftheY

(

1S

)

mesontomuons,themeasuredratioofthecross sec-tionatacenter-of-massof13 TeV tothatat8 TeV is1

.

14

±

0

.

32, wheretheuncertaintyincludesboththestatisticalandsystematic components.The pythia generator predicts aratio of2.1for DPS production,and1.6fortheSPS production.Takingthefractionof the DPSmechanism in the total cross section fDPS

= (

39

±

14

)

%

atacenter-of-mass energyof13 TeV,asmeasuredinSection 5.5, thecrosssectionratiopredictedby pythia is1

.

79

±

0

.

27. Combin-ing the uncertainties in quadrature,the prediction iswithin two standarddeviationsofthemeasurement.

Another unbinned extended maximum likelihood ﬁt is per-formed to extract the number of Y

(

1S

)

Y

(

1S

)

events observed in dataaftertheselection.TheY

(

1S

)

Y

(

1S

)

unweightedsignalyieldis obtainedfrom a fit where all observed events have a weight of 1.0. For this fit, a separate signal shape is determined by fitting the

m

12 andm34 distributionsin theunweighted simulation.The

absenceofweighting doesnotsigniﬁcantly modifythesignal dis-tribution.Theunweighted eventyieldsare givenforallprocesses inTable2.Thereisnoevidenceforthesimultaneousproductionof twoexcitedstatesoftheY meson,butexcesseswithasigniﬁcance

Table 2

Theunweightednumberofeventsforeachofthe pro-cesses from the ﬁt to the m12 and m34 distributions

withoutacceptancenoreﬃciencycorrections. Process Uncorrected yield Y(1S)+Y(1S) 111±16 Y(2S)+Y(2S) 3.6+4.4 −3.6 Y(3S)+Y(3S) 1.1+₋11..41 Y(1S)+combinatorial 166±33 Y(2S)+combinatorial 25±18 Y(3S)+combinatorial 1.1+11 −1.1 Y(2S)+Y(1S) 19±10 Y(3S)+Y(1S) 17±11 Combinatorial+combinatorial 561±41

Fig. 3. Thenumberofdataeventsineach0.6 GeV×0.6 GeV binofthem12vs.m34

distributionisshown.Theresultsofthemaximum-likelihoodﬁttothe signal+back-groundmodelaregivenbythecolors,usingthecolorscaletotherightoftheplot.

lowerthantwostandarddeviationsindicatethepossiblepresence ofY

(

1S

)

Y

(

2S

)

andY

(

1S

)

Y

(

3S

)

events.Thenumberofeventsfrom datainthe

m

12vs.

m

34distributionisshowninFig.3,alongwith

the results ofthe ﬁt to the signal

+

backgroundmodel, using the colorscaletotherightoftheplot.

5.4. Effect of the polarization

Theacceptanceandeﬃciencycorrectionshavebeencomputed assuming negligiblepolarization ofthe Y

(

1S

)

mesons. Adifferent assumptionon thepolarization canchange themeasured ﬁducial cross section.The polarization of theY

(

1S

)

statesaffects the an-gular distributions ofthe leptons produced inthe Y

(

1S

)

→

μ+μ− decaysthroughthefollowingformula [34]:

d2N

d cos

θ

d

φ

∝

1 3

+ λ

θ

(

1

+ λ

θcos2

θ

+ λ

φsin2

θ

cos 2

φ

+ λ

θ φsin 2

θ

cos

φ),

where

θ

and

φ

arethepolarandazimuthalangles,respectively,of the positively charged muon withrespect to the z axis of a po-larization frame,and

λ

θ,

λ

φ,and

λ

θ φ arethe angulardistribution

parameters.Toestimatetheeffectofthepolarizationonthe mea-surement of the Y

(

1S

)

Y

(

1S

)

ﬁducial cross section, we choose to usethehelicityframe,wherethepolaraxiscoincideswiththe di-rectionoftheY

(

1S

)

momentum.Measurementsperformedby the CMSandLHCbCollaborationsonsingleY productionindicate com-patibilityofalltheangulardistributionparameterswithzeroover a large phase space [35,36].However, the same maynot be true for Y

(

1S

)

pair production. To estimate the effect of polarization on the Y

(

1S

)

pairproduction cross section, simulated eventsare

(6)

Table 3

VariationofthemeasuredﬁducialY(1S)pairproductioncrosssectionforseveralλθ

coeﬃcientvalues.

λθ −1.0 −0.5 −0.3 −0.1 +0.1 +0.3 +0.5 +1.0

σﬁd −60% −22% −12% −3.7% +3.4% +9.4% +14% +25%

reweightedtohavetheangulardistributionscorrespondingto var-ious

λ

θ values,withoutchanging theoverall simulatedyield.The

sameeﬃciencyandacceptancecorrectionsasinEq. (2) areusedto calculate

N

corr_for_these_different_polarization_scenarios._The

varia-tionsinthemeasuredY

(

1S

)

pairproductioncrosssectionaregiven fordifferent

λ

θ coeﬃcientsinTable 3. Theeffectofdifferent

po-larizationscanbesubstantial,changingthemeasuredcrosssection by

−

60 to

+

25%.

5.5. Measurement of the DPS-to-inclusive fraction

The DPSandSPS mechanismslead to differentkinematic dis-tributions fortheY

(

1S

)

Y

(

1S

)

events.TheDPSproductionis char-acterizedby a larger separation inrapidity between the mesons,

|

y

(

Y

(

1S

),

Y

(

1S

))

|

,astheyarelargelyuncorrelated,andbyalarger invariantmassofthemesonpairs,

m

Y(1S)Y(1S).Thedistributionsof

φ (

Y

(

1S

),

Y

(

1S

))

,

R

(

Y

(

1S

),

Y

(

1S

))

, and pT

(

Y

(

1S

)

Y

(

1S

))

also

dif-ferfortheSPSandDPSmechanisms,buttheyareverysensitiveto thechoiceofmodelparametersinthesimulationandare subject to large theoretical uncertainties [37]. Measuring the Y

(

1S

)

Y

(

1S

)

ﬁducialcrosssectioninbinsof

|

y

(

Y

(

1S

),

Y

(

1S

))

|

orof

m

Y(1S)Y(1S)

can give a measurement of the fractionof DPSevents, fDPS,

de-ﬁnedas: fDPS

=

σ

DPS ﬁd

σ

SPS ﬁd

+

σ

ﬁdDPS

,

(4) where

σ

DPS ﬁd and

σ

SPS

ﬁd are,respectively,theDPSandSPScross

sec-tionsinthefiducialregion.Wemeasurethefiducialcrosssection infivebins of

|

y

(

Y

(

1S

),

Y

(

1S

))

|

andﬁvebinsof

m

Y(1S)Y(1S).The

signal and background models are the same as forthe inclusive measurement, except that the width of the function describing the Y

(

1S

)

invariant mass shape is allowed to ﬂoat between its best-ﬁtvalues fortheinclusive selectionand fortheselection in therelevantexclusive bin.Thisallowsfora potentialdegradation (improvement)of the muon momentum resolutionat high(low) pseudorapidityto betakenintoaccount,since themuon pseudo-rapidityiscorrelatedwithboth

|

y

(

Y

(

1S

),

Y

(

1S

))

|

and

m

Y(1S)Y(1S).

The systematic uncertainties are identical to those presented in Section5.2.

The extracted ﬁducial cross sections as a function of

|

y

(

Y

(

1S

),

Y

(

1S

))

|

andmY(1S)Y(1S) are compared to the expected

distributionsforSPSandDPSproduction,asobtainedinthe ﬁdu-cialregionusing pythia fortheDPSprocess,andfrom HELAC-Onia with the NLO* CSM predictions for the SPS process. The frac-tion fDPS is measured with a binned maximum-likelihood ﬁt of

thesetwo simulateddistributions withﬂoating normalizationsto themeasured ﬁducialcrosssectionsinbins of

|

y

(

Y

(

1S

),

Y

(

1S

))

|

and

m

Y(1S)Y(1S).Asdeterminedfrompseudo-experiments, thebest

precision is expected to be achieved using

|

y

(

Y

(

1S

),

Y

(

1S

))

|

. Theoreticaluncertainties comingfromthechoiceofparton distri-bution functionsandthefactorizationand renormalizationscales are taken into account for both the SPS and DPS predicted dis-tributions. The fraction fDPS ismeasured to be

(

39

±

14

)

% using

|

y

(

Y

(

1S

),

Y

(

1S

))

|

asthe discriminativedistribution. This results includes both statistical and systematic uncertainties, where the former strongly dominates. The result using mY(1S)Y(1S) is

com-patible with this measurement, but with much lower precision:

(

27

±

22

)

%. The uncertainties are strongly dominated by the

Fig. 4. Measuredﬁducialcrosssection(blackdots)inbinsof|y(Y(1S),Y(1S))| (up-per) ormY(1S)Y(1S) (lower). Thelastbinincludestheoverﬂow.The SPSand DPS

distributionspredictedfromsimulationareoverlaidusingthe fDPSvalueextracted

fromtheﬁttothe|y(Y(1S),Y(1S))|distribution.TheshadedareasaroundtheSPS andDPSpredictionsindicatethetheoreticaluncertainties,whichareoftensmaller thanthe thicknessofthedashedlines.Theshaded areaaroundthetotal distri-butioncorrespondstotheuncertaintyinthemeasurementof fDPS.Thesolidline

showsthesumoftheSPSandDPScontributionswiththebest-ﬁt fDPS.

uncertainties in the measurements of the cross section in the

|

y

(

Y

(

1S

),

Y

(

1S

))

|

and

m

Y(1S)Y(1S)bins,withtheoretical

uncertain-ties in the predictedSPS and DPSdistributions playing a role at thepercentlevel.Themeasureddifferentialﬁducialcrosssections areshowninFig.4,togetherwiththeSPSandDPSpredictions.

6. Searchforresonances 6.1. Methodology

We search for a narrow excess of events above an expected smooth four-muon invariant mass spectrum. Assuming that the resonantstatedecaysintotwomuonsandaY

(

1S

)

mesonthat fur-therdecaystoapairofmuons,thesignal massresolutioncan be improvedbyusingamass-differenceobservable [38]:

m4μ

=

m4μ

−

mμμ

+

mY(1S)

,

(5)

where

m

4μ isthe invariantmass ofthe fourleptons,

m

μμ the

in-variant massassociatedwiththe Y

(

1S

)

candidate,and

m

Y(1S) the

nominalmassoftheY

(

1S

)

particle(9.46 GeV [31]).Thisestimated mass,denotedasm

4μ,hasaresolutionabout50% betterthanthe

four-muoninvariantmass

m

4μ forsignal events.The

m

4μ andm

4μ

(7)

Theresultsareextractedbyperforminganunbinned maximum-likelihood ﬁt to the

m4μ spectrum. The signal and background

components are modeled by several functional forms in the ﬁt, asdescribedinthenextparagraphs.

Thesignal distributions are parameterizedby the sumof two Gaussian functionswiththe same mean.The parameters are ex-tractedforthefourmasspointsavailableinsimulation.Thesignal modeling needs to be interpolated formasses between16.5and 26 GeV and extrapolated to massesup to 27 GeV to search for narrowresonances withanymassbetween16.5and27 GeV.This isdonebyﬁttingwithpolynomialsthedifferentparametersofthe two Gaussian functionsasa function ofthe generatedresonance mass.Thesameprocedureisrepeatedforeverysignal model.The fullwidthathalfmaximumisabout0.2 GeV foraresonancemass of18 GeV.

The background is separated into two components: the Y

(

1S

)

Y

(

1S

)

process,whichwasthesignalinSection5andis char-acterizedbyasharprisingedgeinthem

4μ spectrumattwicethe

Y

(

1S

)

meson mass, and the combinatorial background, which is describedbyasmoothfunctionasexplainedbelow.

Them

4μ spectrumfortheY

(

1S

)

Y

(

1S

)

processisobtainedfrom

simulation,andis modeled asthe product ofa sigmoid function andanexponential functionwithanegative exponent.The nomi-nalmodel for the Y

(

1S

)

Y

(

1S

)

backgroundis takenas an average between the DPS and SPS templates, which is consistent with the measurement of the DPS fraction presented in Section 5.3. Fig. 5 shows the m

4μ models obtained from simulated DPS and

SPS events, together withthe average ﬁt model. The number of Y

(

1S

)

Y

(

1S

)

eventsinthe signalregion isextracted, asdetailedin Section 5, using the selection designed for the resonance search and without applying the acceptance and eﬃciency corrections fromEq. (2).Inthiscase,onlyeventswith13

<

m

4μ

<

28 GeV are

retainedandnorapiditycriteriaareappliedforthereconstructed Y

(

1S

)

candidates.Theyieldismeasuredtobe78

±

13 events.The requirementthatthemassofadimuonpairiscompatiblewiththe massofaY

(

1S

)

mesonwithintwostandarddeviationsisenforced in the resonance search but is not applied to extract the yield becausethe 2D ﬁt relies onthe mass tails toestimate the com-binatorialbackground.Sincetheeﬃciencyofthiscriterion is95% inboth the SPS andDPS Y

(

1S

)

Y

(

1S

)

simulations, the Y

(

1S

)

Y

(

1S

)

yield in the signal region is corrected to 74

±

13. The normal-ization of the Y

(

1S

)

pair production process and its uncertainty are extracted fromthe same data asin the signal region of the resonancesearch,butthisdoesnotlead toa signiﬁcant overcon-straintoftheuncertaintyinthemaximum-likelihoodﬁtofthem

4μ

distributionbecausethe lattercan determinethe Y

(

1S

)

pair nor-malizationonlywithpoorprecision.

Them

4μspectrumforthecombinatorialbackgroundisobtained

intheﬁttothedatainthesignalregion.Severalgenericfunctions areusedtoparameterizethissmoothbackground:

•

Chebychevpolynomialsofvariousorders;

•

thesumofaGaussianfunctionandaChebychevpolynomial;

•

thesumofaBreit–WignerfunctionandaChebychev

polyno-mial.

The widthsof the Gaussian andBreit–Wigner functionsare con-strainedtobeabove2 GeV to avoidfittingnarrowstructuresdue tostatisticalfluctuations.We verify,usingacontrol regionwhere the vertex fit

χ

2 _probability _of _the _four _muons _is _in _the _range

10−10–10−3,thatthesethreefunctionalformsdescribethesmooth

m4μ spectrum of the combinatorial background with a good

χ

2

probability.Muonswithavertexprobabilityinthisrangearelikely tobeassociatedwithprocessesfromthesameprimaryvertex,but canoriginatefromdecaysinﬂightordisplacedsecondaryvertices. ThiscontrolregionisshowninFig.6forillustrativepurposes.The

Fig. 5. Distributionsofm4μforsimulatedY(1S)Y(1S)events.Thedashedlinesare

thebest-ﬁtmodelsfortheSPSandDPSsimulations.Thesolidlineistheaverageof theSPSandDPSmodels,whichistakenasthenominalmodelfortheY(1S)Y(1S) backgroundintheresonancesearch.

Fig. 6. Distributionsofm4μ for thecombinatorialbackgroundinacontrol region

withthevertexﬁt

χ

2_probability_of_the_four_muons_in_the_range₁₀−10_–10−3_._The

parametersobtainedinthisfitarenotusedasaninputforthefitinthesignal region.Thefunctionalformsforthecombinatorialbackgroundshownbythelines areallconsideredaspossibleshapesforthebackgroundmodelinthelikelihoodfit. Theorderofthepolynomialsisindicatedinparenthesesinthelegend.

parameters ofthe functionsdetermined fromtheﬁt arenotused in the signal region, where the parameters of the combinatorial background,aswellasthechoiceofthefunctionalform,arefreely ﬂoating.

6.2. Systematic uncertainties

The systematic uncertainties are to a large extent similar to those used in the measurement of the Y

(

1S

)

Y

(

1S

)

cross section andintroduced inSection 5.2.Inthissection,onlythedifferences arehighlighted.Theyarisefromslightlydifferentselectioncriteria, a differentchoice of observable, thetreatment of the Y

(

1S

)

Y

(

1S

)

processasabackground,andtheintroductionofanewsignal pro-cess.

Theuncertaintypermuoninthemuonidentiﬁcationand track-ing is increased from 0.5% to 1% because poorly reconstructed muonswith pT

<

3

.

5 GeV inthe barrelare includedinthe

reso-nancesearchtoincreasethesignalacceptanceforlightresonances. Inaddition,intheresonancesearch,thesignalisaffectedbya1% yield uncertainty related to the requirementthat the Y

(

1S

)

(8)

can-Fig. 7. Them4μ distributionfromdataandtheresultsoftheﬁtintheresonance

search.Anexamplesignalisshownforthetetraquarkmodelwithamassof19 GeV, whichhasasigniﬁcanceofaboutonestandarddeviation.

didatehasan invariant masscompatiblewith thenominal Y

(

1S

)

meson mass within two standard deviations. This uncertainty is determined by comparing the dimuon invariant mass resolution distributions in Y

(

1S

)

Y

(

1S

)

simulated eventsand inY

(

1S

)

events indata.Themodelingofthesignalprocesswitharesonancemass other than those for which simulated samples were generated leadstoa2%uncertaintyinthesignalnormalizationforthe reso-nancesearch.

The discreteprofiling method [39] isusedto modelthe com-binatorial background.This allows the choice ofthe fit functions amongthoseprovidedtobeconsideredasadiscretenuisance pa-rameter.Theparametersofthesefitfunctionsarefreelyfloating.

The normalization of the Y

(

1S

)

Y

(

1S

)

background in the reso-nance search is extracted from the 2D unbinned ﬁt to the in-variant massof the dimuon pairs in the Y

(

1S

)

massregion. The uncertainty inthe yield obtained fromthe ﬁt isconsidered asa log-normaluncertaintyintheﬁttothe

m4μ distribution.The

m4μ

distributionof theY

(

1S

)

Y

(

1S

)

background isallowed to ﬂoat be-tweenthepredictionsfortheSPSandDPSsimulations.

Uncertaintiesinthem

4μdistributionoftheresonantsignaltake

into account the limited size of the simulated samples, andthe limited precision of the description of the signal formasses not availableinsimulations.Theuncertaintyinthemeanmassofthe signalis 0.2%,correspondingto theuncertaintyinthemuon mo-mentum scale.The other parameters describingthe shape ofthe signal havean uncertainty between5 and15%, whichleads to a combinedimpactontheﬁnalupperlimitsoflessthan2%.

Theuncertaintyinthe Y

(

1S

)

dimuon branchingfractionisnot considered,sincethelimitsaresetontheproductoftheresonance productioncrosssectionandits branchingfractiontofourmuons viaanintermediateY

(

1S

)

resonance.

6.3. Results

The binned

m4μ distribution in the signal region of the

res-onance search is shown in Fig. 7. The background andexample signalcomponentsareshownusingtheir best-ﬁtshapesand nor-malizations. Using the number ofY

(

1S

)

Y

(

1S

)

events observed in dataasareference, aresonancewitha massaround19 GeV and having a similar production cross section andbranching fraction tofourmuonsastheY

(

1S

)

Y

(

1S

)

production,wouldproduceabout 100eventsinoursample, giventhesimilarity betweenthe kine-matic distributions of both processes. No signiﬁcant narrow ex-cessofeventsisobservedabove thebackgroundexpectation.The largestexcess isobserved fora resonancemass of25.1 GeV, and

has a local signiﬁcance of 2.4 standard deviations for the scalar signalhypothesis.

Upperlimitsontheproductoftheproductioncrosssectionofa resonanceandthebranchingfractiontoaﬁnalstateoffourmuons viaanintermediateY

(

1S

)

resonancearesetat95%conﬁdencelevel (CL)usingthemodiﬁedfrequentistconstructionCLsinthe

asymp-toticapproximation [40–44],separatelyforeachsignalmodel.The upperlimitsare extractedusingunbinneddistributions.Thecross section is deﬁned in the entire phase space without ﬁducial re-quirements,andthebranching fractionusedistheproductofthe branching fraction of the resonant state to a Y

(

1S

)

meson and two muons, and the branching fraction of the Y

(

1S

)

meson to twomuons. Massesbetween17.5and19.5 GeV areprobedinthe context of the tetraquark search, using the bottomonium model, whereasthelimitsintheextendedmassrange16.5–27 GeV areset forthegeneric search,usingthe JHUGen models.The correspond-ing upper limits are given in Fig. 8. They range between 5 and 380 fb,depending on themassandsignal model.The patternsin thelimitsarebroaderforthespin-2signalthanforthescalarand pseudoscalar states because the signal is characterized by softer andmoreforwardmuons,leadingtoaworse

m4μ resolution. 7. Summary

The crosssection forY

(

1S

)

pairproductionismeasuredinthe ﬁducial regionwhereboth Y

(

1S

)

mesonshavean absolute rapid-itybelow2.0.Themeasurementisperformedusingproton-proton collision data collected at a center-of-mass energy of 13 TeV by the CMSdetectorin2016andcorresponding to anintegrated lu-minosity of 35.9 fb−1. Assuming that the Y

(

1S

)

mesons are pro-ducedunpolarized,theﬁducialY

(

1S

)

pairproductioncrosssection is determinedto be 79

±

11(stat)

±

6(syst)

±

3

(B)

pb, wherethe last uncertaintycomesfromtheuncertaintyin theY

(

1S

)

dimuon branching fraction.TheresultcanchangeiftheY

(

1S

)

mesonsare produced witha nonzero polarization. Changingthe polarization coeﬃcient

λ

θ from

−

1 to

+

1,theresultingY

(

1S

)

pairproduction

crosssectionmeasurementvariesby

−

60 to

+

25%.

Thecontributionofdouble-partonscatteringtothetotal inclu-siveY

(

1S

)

pairproductioncrosssectionisdeterminedfortheﬁrst time. Itismeasured tobe (39

±

14)%inthesame ﬁducialregion asdescribedabove,wheretheuncertaintyincludesbothstatistical andsystematiccomponents,withthestatisticaluncertainty domi-nating.

Theresultsofasearcharealsopresentedforalightnarrow res-onance,suchasatetraquarkoraboundstatebeyond-the-standard model,decayingtoaY

(

1S

)

andapairofopposite-signmuons.No excessofeventscompatiblewitha signalisobservedinthe four-muon invariant mass spectrum. Upper limits at 95% conﬁdence levelontheproductofthesignalcrosssectionandbranching frac-tiontofourmuonsviaanintermediateY

(

1S

)

resonancearesetfor differentsignalmodels,expandingthekinematicandmass cover-ageofprevioussearches.

Declarationofcompetinginterest

Theauthorsdeclarethattheyhavenoknowncompeting ﬁnan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto inﬂuencetheworkreportedinthispaper.

Acknowledgements

WecongratulateourcolleaguesintheCERNaccelerator depart-ments for the excellent performance of the LHC and thank the technical andadministrativestaffs atCERNand atother CMS in-stitutes for their contributions to the success of the CMS effort. Inaddition,wegratefullyacknowledgethecomputingcentersand

(9)

Fig. 8. Upperlimitsat95%CL ontheproductofthecrosssectionandbranchingfractionforatetraquark(upperleft),scalar(upperright),pseudoscalar(lowerleft),and spin-2(lowerright)states.Thesymbol

σ

denotestheproductioncrosssectionoftheresonance,andthesymbolBdenotestheproductofthebranchingfractionforthe decayoftheresonancetoaY(1S)mesonandtwomuons,andtheY(1S)mesondimuonbranchingfraction.Thelinewiththepointsonitshowstheobservedupperlimits andthethinredlineisthemedianoftheexpectedupperlimits.Theinner(green)bandandtheouter(yellow)bandindicatetheregionscontaining68and95%,respectively, ofthedistributionoflimitsexpectedunderthebackground-onlyhypothesis.

personneloftheWorldwideLHCComputingGridfordeliveringso effectivelythecomputinginfrastructure essential toour analyses. Finally, we acknowledge the enduring support for the construc-tionandoperationofthe LHCandtheCMSdetectorprovided by thefollowingfundingagencies: BMBWFandFWF(Austria);FNRS andFWO (Belgium); CNPq, CAPES, FAPERJ,FAPERGS, andFAPESP (Brazil); MES (Bulgaria); CERN; CAS, MOST, and NSFC (China); COLCIENCIAS (Colombia); MSES and CSF (Croatia); RPF (Cyprus); SENESCYT (Ecuador); MoER, ERC IUT, PUT and ERDF (Estonia); AcademyofFinland,MEC,andHIP(Finland);CEAandCNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany); GSRT (Greece); NK-FIA (Hungary); DAE and DST (India); IPM (Iran); SFI (Ireland); INFN(Italy);MSIPandNRF(RepublicofKorea);MES(Latvia);LAS (Lithuania);MOEandUM(Malaysia); BUAP,CINVESTAV,CONACYT, LNS,SEP,andUASLP-FAI(Mexico);MOS(Montenegro);MBIE(New Zealand); PAEC (Pakistan); MSHE and NSC (Poland); FCT (Portu-gal);JINR(Dubna);MON,ROSATOM, RAS,RFBR,andNRCKI (Rus-sia);MESTD(Serbia);SEIDI,CPAN,PCTI,andFEDER(Spain);MoSTR (Sri Lanka); Swiss Funding Agencies (Switzerland); MST (Taipei); ThEPCenter,IPST,STAR,andNSTDA(Thailand);TÜBITAK andTAEK

(Turkey); NASU (Ukraine); STFC (United Kingdom);DOE andNSF (USA).

Individuals have received support from the Marie-Curie pro-gramandtheEuropeanResearchCouncilandHorizon2020Grant, contractNos.675440,752730, and765710 (EuropeanUnion);the Leventis Foundation; the AlfredP. Sloan Foundation;the Alexan-der vonHumboldt Foundation;theBelgianFederalScience Policy Oﬃce; the Fonds pour la Formation à la Recherche dans l’In-dustrieetdans l’Agriculture(FRIA-Belgium); the Agentschapvoor Innovatie door Wetenschap en Technologie (IWT-Belgium); the F.R.S. - FNRS and FWO (Belgium) under the “Excellence of Sci-ence – EOS” – be.h project n. 30820817; the Beijing Municipal ScienceandTechnology Commission,No.Z191100007219010;the Ministry of Education, Youth and Sports (MEYS) of the Czech Republic;theDeutscheForschungsgemeinschaft(DFG)under Ger-many’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306; the Lendület (“Momentum”) Program and the János Bolyai Research Scholarship of the Hungarian Academy of Sci-ences, the New National Excellence Program ÚNKP, the NK-FIA research grants 123842, 123959, 124845, 124850, 125105,

(10)

128713, 128786, and 129058 (Hungary); the Council of Science andIndustrialResearch,India; theHOMING PLUSprogramofthe Foundation for Polish Science, cofinanced from European Union, Regional Development Fund, the Mobility Plus program of the Ministry of Science and Higher Education, the National Science Center (Poland), contractsHarmonia 2014/14/M/ST2/00428, Opus 2014/13/B/ST2/02543, 2014/15/B/ST2/03998, and 2015/19/B/ST2/ 02861,Sonata-bis2012/07/E/ST2/01406;theNationalPriorities Re-search Program by Qatar National Research Fund; the Ministry of Science and Education, grant no. 14.W03.31.0026(Russia); the ProgramaEstatalde Fomentode laInvestigación Científicay Téc-nica de Excelencia María de Maeztu, grant MDM-2015-0509 and theProgramaSeveroOchoadelPrincipadode Asturias;theThalis andAristeiaprograms cofinancedbyEU-ESF andtheGreek NSRF; theRachadapisekSompotFundforPostdoctoralFellowship, Chula-longkornUniversityandtheChulalongkornAcademicintoIts 2nd Century Project Advancement Project (Thailand);The Kavli Foun-dation;the Nvidia Corporation; the SuperMicroCorporation; The WelchFoundation,contractC-1845;andtheWestonHavens Foun-dation(USA).

References

[1]J.-P.Lansberg,New observablesininclusiveproduction ofquarkonia,arXiv: 1903.09185,2019.

[2] CMSCollaboration,ObservationofY(1S)pairproductioninproton-proton col-lisionsat√s=8 TeV,J.HighEnergyPhys.05(2017)013,https://doi.org/10. 1007/JHEP05(2017)013,arXiv:1610.07095.

[3] Y.Chen,S.Wu,Productionoffour-quarkstateswithdoubleheavyquarksat LHC,Phys.Lett.B705(2011)93,https://doi.org/10.1016/j.physletb.2011.09.096, arXiv:1101.4568.

[4] A.V. Berezhnoy, A.V. Luchinsky, A.A. Novoselov, Tetraquarks composed of 4 heavy quarks, Phys. Rev. D 86 (2012) 034004, https://doi.org/10.1103/ PhysRevD.86.034004,arXiv:1111.1867.

[5] B.A. Dobrescu, C. Frugiuele, Hidden GeV-scaleinteractions of quarks,Phys. Rev.Lett.113(2014)061801,https://doi.org/10.1103/PhysRevLett.113.061801, arXiv:1404.3947.

[6] M. Karliner,S. Nussinov,J.L. Rosner,QQQ¯Q states:¯ masses,production,and decays,Phys. Rev.D95(2017) 034011, https://doi.org/10.1103/PhysRevD.95. 034011,arXiv:1611.00348.

[7] W. Chen,H.-X.Chen,X.Liu,T.G.Steele,S.-L.Zhu,Huntingforexoticdoubly hidden-charm/bottomtetraquarkstates,Phys.Lett.B773(2017)247,https:// doi.org/10.1016/j.physletb.2017.08.034,arXiv:1605.01647.

[8] Z.-G. Wang,Analysis ofthe QQQ¯Q tetraquark¯ states with QCD sum rules, Eur.Phys.J.C77(2017)432,https://doi.org/10.1140/epjc/s10052-017-4997-0, arXiv:1701.04285.

[9] J.Wu,Y.-R.Liu,K.Chen,X.Liu,S.-L.Zhu,Heavy-ﬂavoredtetraquarkstateswith the QQQ¯Q conﬁguration,¯ Phys. Rev.D97(2018) 094015,https://doi.org/10. 1103/PhysRevD.97.094015,arXiv:1605.01134.

[10] M.N.Anwar,J.Ferretti,F.-K.Guo,E.Santopinto,B.-S.Zou,Spectroscopyand decays ofthefully-heavytetraquarks,Eur.Phys.J.C78(2018)647,https:// doi.org/10.1140/epjc/s10052-018-6073-9,arXiv:1710.02540.

[11] C. Hughes, E. Eichten, C.T.H. Davies, Searching for beauty-fully bound tetraquarksusinglatticenonrelativisticQCD,Phys.Rev.D97(2018)054505, https://doi.org/10.1103/PhysRevD.97.054505,arXiv:1710.03236.

[12] A.Esposito,A.D.Polosa,Abbb¯b di-bottomonium¯ attheLHC?,Eur.Phys.J.C78 (2018)782,https://doi.org/10.1140/epjc/s10052-018-6269-z,arXiv:1807.06040. [13] Y. Bai, S. Lu, J. Osborne, Beauty-full tetraquarks, Phys. Lett. B 798(2019)

134930,https://doi.org/10.1016/j.physletb.2019.134930,arXiv:1612.00012. [14] LHCb Collaboration, Search for beautiful tetraquarks in the Y(1S)μ+μ−

invariant-massspectrum,J.HighEnergyPhys.10(2018)086,https://doi.org/ 10.1007/JHEP10(2018)086,arXiv:1806.09707.

[15] CMS Collaboration, The CMStrigger system, J. Instrum. 12 (2017)P01020, https://doi.org/10.1088/1748-0221/12/01/P01020,arXiv:1609.02366. [16] CMSCollaboration,TheCMSexperimentattheCERNLHC,J.Instrum.3(2008)

S08004,https://doi.org/10.1088/1748-0221/3/08/S08004.

[17] CMSCollaboration,PerformanceoftheCMSmuondetectorandmuon recon-structionwithproton-protoncollisionsat√s=13TeV,J.Instrum.13(2018) P06015,https://doi.org/10.1088/1748-0221/13/06/P06015,arXiv:1804.04528. [18] T.Sjöstrand, S.Ask,J.R.Christiansen,R.Corke, N.Desai,P.Ilten,S.Mrenna,

S.Prestel,C.O.Rasmussen,P.Z.Skands,AnintroductiontoPYTHIA8.2, Com-put.Phys.Commun.191(2015)159,https://doi.org/10.1016/j.cpc.2015.01.024, arXiv:1410.3012.

[19] R. Baier,R. Ruckl, Hadronic production of J/ and ϒ: transverse momen-tumdistributions,Phys.Lett.B102(1981)364,https://doi.org/10.1016/0370 -2693(81)90636-5.

[20] J.-P.Lansberg,H.-S.Shao,ProductionofJ/ψ+ηcversusJ/ψ+J/ψattheLHC:

importanceofreal

α

5

s corrections,Phys.Rev.Lett.111(2013)122001,https://

doi.org/10.1103/PhysRevLett.111.122001,arXiv:1308.0474.

[21] J.-P.Lansberg,H.-S.Shao, J/ψ-pairproductionatlargemomenta:indications for doubleparton scatteringsand largeα5

s contributions, Phys. Lett.B751

(2015)479,https://doi.org/10.1016/j.physletb.2015.10.083,arXiv:1410.8822. [22] H.-S. Shao, HELAC-Onia: an automatic matrix element generator for heavy

quarkoniumphysics,Comput.Phys.Commun.184(2013)2562,https://doi.org/ 10.1016/j.cpc.2013.05.023,arXiv:1212.5293.

[23] H.-S.Shao,HELAC-Onia 2.0:anupgradedmatrix-elementandevent genera-torfor heavyquarkoniumphysics,Comput.Phys.Commun.198(2016)238, https://doi.org/10.1016/j.cpc.2015.09.011,arXiv:1507.03435.

[24] Y.Gao,A.V.Gritsan,Z.Guo,K.Melnikov,M.Schulze,N.V.Tran,Spin determina-tionofsingle-producedresonancesathadroncolliders,Phys.Rev.D81(2010) 075022,https://doi.org/10.1103/PhysRevD.81.075022,arXiv:1001.3396. [25] S.Bolognesi,Y.Gao,A.V.Gritsan,K.Melnikov,M.Schulze,N.V.Tran,A.

Whit-beck,On the spin and parity of a single-produced resonance at the LHC, Phys. Rev.D86 (2012) 095031, https://doi.org/10.1103/PhysRevD.86.095031, arXiv:1208.4018.

[26] I.Anderson,S.Bolognesi,F.Caola,Y.Gao,A.V.Gritsan,C.B.Martin,K.Melnikov, M.Schulze,N.V.Tran,A.Whitbeck,Y.Zhou,ConstraininganomalousHVV inter-actionsatprotonandleptoncolliders,Phys.Rev.D89(2014)035007,https:// doi.org/10.1103/PhysRevD.89.035007,arXiv:1309.4819.

[27] A.V.Gritsan,R.Röntsch,M.Schulze,M.Xiao,ConstraininganomalousHiggs bosoncouplingstotheheavyﬂavorfermionsusingmatrixelementtechniques, Phys. Rev.D94(2016) 055023, https://doi.org/10.1103/PhysRevD.94.055023, arXiv:1606.03107.

[28] CMSCollaboration,Eventgeneratortunesobtainedfromunderlyingeventand multipartonscatteringmeasurements,Eur.Phys.J.C76(2016)155,https:// doi.org/10.1140/epjc/s10052-016-3988-x,arXiv:1512.00815.

[29] S.Agostinelli,et al.,GEANT4,GEANT4—asimulationtoolkit,Nucl.Instrum. MethodsA506(2003)250,https://doi.org/10.1016/S0168-9002(03)01368-8. [30] CMS Collaboration,Particle-ﬂowreconstruction andglobalevent description

withtheCMSdetector,J.Instrum.12(2017)P10003,https://doi.org/10.1088/ 1748-0221/12/10/P10003,arXiv:1706.04965.

[31] ParticleDataGroup,M.Tanabashi,etal.,Reviewofparticlephysics,Phys.Rev. D98(2018)030001,https://doi.org/10.1103/PhysRevD.98.030001.

[32] R.A.Fisher,Ontheinterpretationof

χ

2_from_contingency_tables,_and_the

cal-culationofp,J.R.Stat.Soc.85(1922)87,https://doi.org/10.2307/2340521. [33] CMSCollaboration,CMSLuminosityMeasurementsforthe2016DataTaking

Period,2017,CMSPhysicsAnalysisSummary,CMS-PAS-LUM-17-001,https:// cds.cern.ch/record/2257069.

[34] P.Faccioli,C.Lourenco,J.Seixas,H.K.Wohri,Towardstheexperimental clar-iﬁcationofquarkoniumpolarization, Eur.Phys. J. C69(2010) 657,https:// doi.org/10.1140/epjc/s10052-010-1420-5,arXiv:1006.2738.

[35] LHCb_√ Collaboration,Measurementoftheϒpolarizationsin pp collisionsat s=7 and8TeV,J.HighEnergyPhys.12(2017)110,https://doi.org/10.1007/ JHEP12(2017)110,arXiv:1709.01301.

[36] CMSCollaboration,MeasurementoftheY(1S),Y(2S)andY(3S)polarizations inpp collisions at√s=7 TeV,Phys. Rev.Lett.110(2013)081802,https:// doi.org/10.1103/PhysRevLett.110.081802,arXiv:1209.2922.

[37] A.D.Cook,Double J/ψproductionatLHCb,PhDthesis,2015,https://cds.cern. ch/record/2140033.

[38] CMSCollaboration,SearchforHiggsbosonpairproductioninthe

γ γ

bb ﬁnal stateinppcollisionsat√s=13TeV,Phys. Lett.B788(2019)7,https:// doi.org/10.1016/j.physletb.2018.10.056,arXiv:1806.00408.

[39] P.D.Dauncey,M.Kenzie,N.Wardle,G.J.Davies,Handlinguncertaintiesin back-groundshapes, J. Instrum. 10 (2015) P04015, https://doi.org/10.1088/1748 -0221/10/04/P04015,arXiv:1408.6865.

[40] TheATLASCollaboration,TheCMSCollaboration,TheLHCHiggsCombination Group,Procedure for the LHCHiggsboson searchcombination inSummer 2011, Technical Report CMS-NOTE-2011-005. ATL-PHYS-PUB-2011-11, 2011, https://cds.cern.ch/record/1379837.

[41] S.Chatrchyan,etal.,CMS,Combinedresultsofsearchesforthestandardmodel Higgsbosoninppcollisionsat√s=7 TeV,Phys.Lett.B710(2012)26,https:// doi.org/10.1016/j.physletb.2012.02.064,arXiv:1202.1488.

[42] T.Junk,Conﬁdencelevelcomputationforcombiningsearcheswithsmall statis-tics,Nucl.Instrum.MethodsA434(1999)435,https://doi.org/10.1016/S0168 -9002(99)00498-2,arXiv:hep-ex/9902006.

[43] A.L.Read,Presentationofsearchresults:theC Lstechnique,J.Phys.G28(2002)

2693,https://doi.org/10.1088/0954-3899/28/10/313.

[44] G.Cowan,K.Cranmer,E.Gross,O.Vitells,Asymptoticformulaefor likelihood-based tests of new physics, Eur. Phys. J. C 71 (2011) 1554, https://doi. org/10.1140/epjc/s10052-011-1554-0,arXiv:1007.1727,https://doi.org/10.1140/ epjc/s10052-013-2501-z.

(11)

TheCMSCollaboration

A.M. Sirunyan

†

,

A. Tumasyan

YerevanPhysicsInstitute,Yerevan,Armenia

W. Adam,

F. Ambrogi,

T. Bergauer,

M. Dragicevic,

J. Erö,

A. Escalante Del Valle,

M. Flechl,

R. Frühwirth

1

,

M. Jeitler

1

,

N. Krammer,

I. Krätschmer,

D. Liko,

T. Madlener,

I. Mikulec,

N. Rad,

J. Schieck

1

,

R. Schöfbeck,

M. Spanring,

W. Waltenberger,

C.-E. Wulz

1

,

M. Zarucki

InstitutfürHochenergiephysik,Wien,Austria

V. Drugakov,

V. Mossolov,

J. Suarez Gonzalez

InstituteforNuclearProblems,Minsk,Belarus

M.R. Darwish,

E.A. De Wolf,

D. Di Croce,

X. Janssen,

A. Lelek,

M. Pieters,

H. Rejeb Sfar,

H. Van Haevermaet,

P. Van Mechelen,

S. Van Putte,

N. Van Remortel

UniversiteitAntwerpen,Antwerpen,Belgium

F. Blekman,

E.S. Bols,

S.S. Chhibra,

J. D’Hondt,

J. De Clercq,

D. Lontkovskyi,

S. Lowette,

I. Marchesini,

S. Moortgat,

Q. Python,

S. Tavernier,

W. Van Doninck,

P. Van Mulders

VrijeUniversiteitBrussel,Brussel,Belgium

D. Beghin,

B. Bilin,

B. Clerbaux,

G. De Lentdecker,

H. Delannoy,

B. Dorney,

L. Favart,

A. Grebenyuk,

A.K. Kalsi,

L. Moureaux,

A. Popov,

N. Postiau,

E. Starling,

L. Thomas,

C. Vander Velde,

P. Vanlaer,

D. Vannerom

UniversitéLibredeBruxelles,Bruxelles,Belgium

T. Cornelis,

D. Dobur,

I. Khvastunov

2

,

M. Niedziela,

C. Roskas,

K. Skovpen,

M. Tytgat,

W. Verbeke,

B. Vermassen,

M. Vit

GhentUniversity,Ghent,Belgium

O. Bondu,

G. Bruno,

C. Caputo,

P. David,

C. Delaere,

M. Delcourt,

A. Giammanco,

V. Lemaitre,

J. Prisciandaro,

A. Saggio,

M. Vidal Marono,

P. Vischia,

J. Zobec

UniversitéCatholiquedeLouvain,Louvain-la-Neuve,Belgium

G.A. Alves,

G. Correia Silva,

C. Hensel,

A. Moraes

CentroBrasileirodePesquisasFisicas,RiodeJaneiro,Brazil

E. Belchior Batista Das Chagas,

W. Carvalho,

J. Chinellato

3

,

E. Coelho,

E.M. Da Costa,

G.G. Da Silveira

4

,

D. De Jesus Damiao,

C. De Oliveira Martins,

S. Fonseca De Souza,

H. Malbouisson,

J. Martins

5

,

D. Matos Figueiredo,

M. Medina Jaime

6

,

M. Melo De Almeida,

C. Mora Herrera,

L. Mundim,

H. Nogima,

W.L. Prado Da Silva,

P. Rebello Teles,

L.J. Sanchez Rosas,

A. Santoro,

A. Sznajder,

M. Thiel,

E.J. Tonelli Manganote

3

,

F. Torres Da Silva De Araujo,

A. Vilela Pereira

UniversidadedoEstadodoRiodeJaneiro,RiodeJaneiro,Brazil

C.A. Bernardes

a

,

L. Calligaris

a

,

T.R. Fernandez Perez Tomei

a

,

E.M. Gregores

b

,

D.S. Lemos,

P.G. Mercadante

b

,

S.F. Novaes

a

,

Sandra

S. Padula

a

a_Universidade_Estadual_Paulista,_São_Paulo,_Brazil b_Universidade_Federal_do_ABC,_São_Paulo,_Brazil

A. Aleksandrov,

G. Antchev,

R. Hadjiiska,

P. Iaydjiev,

M. Misheva,

M. Rodozov,

M. Shopova,

G. Sultanov

InstituteforNuclearResearchandNuclearEnergy,BulgarianAcademyofSciences,Soﬁa,Bulgaria