Evidence for a new structure in the $J/\psi p$ and $J/\psi \bar{p}$ systems in $B_s^0 \to J/\psi p \bar{p}$ decays

An amplitude analysis of flavour-untagged $B_s^0 \to J/\psi p \bar{p}$ decays is performed using a sample of $797\pm31$ decays reconstructed with the LHCb detector. The data, collected in proton-proton collisions between 2011 and 2018, correspond to an integrated luminosity of 9 $\text{fb}^{-1}$. Evidence for a new structure in the $J/\psi p$ and $J/\psi \bar{p}$ systems with a mass of $4337 \ ^{+7}_{-4} \ ^{+2}_{-2}~\text{MeV}$ and a width of $29 \ ^{+26}_{-12} \ ^{+14}_{-14}~\text{MeV}$ is found, where the first uncertainty is statistical and the second systematic, with a significance in the range of 3.1 to 3.7 $\sigma$, depending on the assigned $J^P$ hypothesis.

The observation of pentaquark candidates (P c ) in J/ψp final states produced in Λ 0 b → J/ψpK − decays 1 [1,2] by the LHCb experiment has stimulated interest in exotic spectroscopy. Recently, evidence for a structure in the J/ψΛ invariant-mass spectrum, consistent with a charmonium-like pentaquark with strangeness, was found in Ξ − b → J/ψΛK − decays [3]. The mass of these states is just below threshold for the joint production of a charm baryon and a charm meson, i.e. the Σ c D * and the Ξ c D * thresholds for the J/ψp and the J/ψΛ resonances, respectively. The mass separation from these thresholds might provide useful information for the phenomenological interpretation for these states. Proposed interpretation can be grouped into three classes: QCD-inspired models [4,5], residual hadron-hadron interaction models [6] and rescattering effects particle [7]. Additional measurements in different productions and decay channels are crucial to disentangle the various models [8].
The B 0 s → J/ψpp decay was observed for the first time by the LHCb experiment in 2019 [9]. This channel may have sensitivity to the resonant P c structures [1,2] within the J/ψp invariant-mass range of [4034,4429] MeV. Additionally, it could proceed via an intermediate glueball candidate f J (2220) decaying to pp [10]. Unlike Λ 0 b → J/ψpK − decays receiving a relatively large contribution from the intermediate excited Λ resonances, no conventional states are expected to be produced in the B 0 s decay, offering a clean environment to search for new resonant structures. Baryonic B 0 (s) decays also allow for a study of the dynamics of the baryon-antibaryon system and its characteristic threshold enhancement, the origin of which is still to be understood [11].
In this Letter, an amplitude analysis of B 0 s → J/ψpp decay is presented, including a search for pentaquark and glueball states, using proton-proton (pp) collision data at centre-of-mass energies of 7 TeV, 8 TeV and 13 TeV, corresponding to a luminosity of 9 fb −1 , collected between 2011 and 2018. The measurement is performed untagged, such that decays of B 0 s and B 0 s are not distinguished and analysed together. The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, described in detail in Refs. [12][13][14][15]. The online event selection is performed by a trigger [16], comprising a hardware stage based on information from the muon system which selects J/ψ → µ + µ − decays, followed by a software stage that applies a full event reconstruction. The software trigger relies on identifying J/ψ decays into muon pairs consistent with originating from a B meson decay vertex detached from the primary pp collision point.
Samples of simulated events are used to study the properties of the signal and control channels. The pp collisions are generated using Pythia [17] with a specific LHCb configuration [18]. Decays of hadronic particles and interactions with the detector material are described by EvtGen [19], using Photos [20], and by the Geant4 toolkit [21,22], respectively. The signal B 0 s → J/ψpp decays are generated from a uniform phase space distribution, while the B 0 s → J/ψφ(→ K + K − ) control mode is generated according to the model of Ref. [23].
The event selection follows the same strategy as Ref. [9]. Signal B 0 s candidates are formed from two pairs of oppositely charged tracks. The first pair is required to be consistent with muons originating from a J/ψ meson with a decay vertex significantly displaced from its associated primary pp vertex (PV). For a given particle, the associated PV is the one with the smallest impact parameter χ 2 IP , defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the track under consideration. The second pair is required to be consistent with protons originating from the muon-pair vertex. A kinematic fit [24] to the B 0 s candidate is performed, with the dimuon mass constrained to the known J/ψ mass [25]. The selection is optimised using multivariate techniques [26] trained with simulation and data. Simulated events are weighted such that the distributions of momentum, p, transverse momentum, p T , and number of tracks per event for B 0 s candidates match the B 0 s → J/ψφ control-mode distributions in data. In simulation the particle identification (PID) variables for each charged track are resampled as a function of its p, p T and the number of tracks in the event using Λ + c → pK − π + and D * + → D 0 (→ K − π + )π + calibration samples from data [27]. The selection consists of two boosted decision tree (BDT) classifiers. The first classifier, BDT sel , is a selection trained on B 0 s → J/ψφ simulation and sideband data with the J/ψpp invariant mass above 5450 MeV using the p, p T , and χ 2 IP variables of the B 0 s candidate, the χ 2 probability from the kinematic fit of the candidate, and the impact parameter distances of the two muons. The second classifier, BDT PID , is trained on B 0 s → J/ψpp simulation and sideband data using proton identification variables: the hadron PID from the ring-imaging Cherenkov detectors, the p, p T and χ 2 IP of the protons. The BDT PID output selection criterion is chosen by maximising the figure of merit S 2 /(S + B) 3/2 , where S and B are the signal and background yields in a region of ±10 MeV around the B 0 s mass peak. These are determined from a fit to the J/ψpp invariant-mass distribution in data after the BDT sel selection, multiplied by the efficiency of the BDT PID output requirement, obtained from simulation and from sideband data, respectively.
After applying these selection criteria, a maximum-likelihood fit is performed to the J/ψpp invariant-mass distribution, shown in Fig. 1, yielding 797 ± 31 B 0 s signal decays. The B 0 s signal shape is modelled as the sum of two Crystal Ball [28] functions sharing a common peak position, with asymmetric tails describing radiative and misreconstruction effects. The signal-model parameters are determined from simulation and only the B 0 s peak position is allowed to vary in the fit to data. The combinatorial background is modelled by a first-order polynomial with parameters determined from the fit to data. The B 0 → J/ψpp component has the same shape as the B 0 s signal. The combinatorial-background fraction in the B 0 s signal window of 3σ around the mass peak ([5357, 5378] MeV) is estimated to be (14.9 ± 0.6)%, where σ ≈ 3.5 MeV is the resolution of the reconstructed invariant mass. The m(J/ψp) and m(J/ψp) invariant mass distributions of the reconstructed B 0 s candidates in the B 0 s signal region are shown in the bottom row of Fig. 2 (black dots), where hints of structure in the region around (4.3 -4.4) GeV are present. This Letter investigates the nature of these enhancements, which are not compatible with the pure phase-space hypothesis.
An amplitude analysis of the B 0 s candidates is performed under the assumption of CP symmetry conservation, i.e. the dynamics is the same in B 0 s and B 0 s decays. Three interfering decay sequences are considered in the amplitude model: all followed by a J/ψ → µ + µ − decay. These sequences are labelled as the X, P + c and P − c chains, respectively. Since the data sample is not flavour tagged, the distribution of the candidates in the phase space is by construction symmetric for J/ψp and J/ψp final states, and therefore the analysis is sensitive to the sum of possible contributions from P + c and P − c pentaquark candidates, denoted as P c in the following. Due to the small sample size and since the B 0 s or B 0 s flavour is not identified, there is no sensitivity to different couplings for the P + c and P − c states, which are constrained to be equal, up to a phase difference. The amplitude model is based on the helicity formalism of Refs. [29,30], which defines a consistent framework for propagating spin correlations through relativistic decay chains. To align the spin of the different decay chains, the prescription in Ref. [31] is followed. Details about the amplitude definition are given in the Supplemental material.
Candidates in the B 0 s signal region are used to perform an amplitude fit in the fourdimensional phase space, (m pp , Ω). This phase space is defined by the invariant mass m pp of the pp pair and Ω = (θ p , θ µ , ϕ), where θ p , θ µ are the two helicity angles of the p and the µ − in the X and J/ψ rest frame, respectively, and ϕ is the azimuthal angle between the decay planes, of the µ − µ + and the pp pairs. The distributions of (m pp , cos θ µ , cos θ p , ϕ), together with the m(J/ψp) and m(J/ψp) invariant mass projections, are shown in Fig. 2 for selected candidates.
The amplitude fit minimises the negative log-likelihood function where the total probability density function (PDF) calculated for i th candidate has a signal, P sig , and a background, P bkg , component where β is the fraction of background events observed within the B 0 s signal window. The signal PDF is proportional to the matrix element squared, |M (m pp,i , Ω i | ω) | 2 , and depends on the fit parameters, ω, i.e. the couplings, the masses and the widths, which define the contributing resonances The phase-space element is Φ (m pp,i ) = | p|| q|, where p is the momentum of the X system in the B 0 s rest frame and q is the proton momentum in the X rest frame. The efficiency, (m pp,i , Ω i ), is included in the PDF, and is parameterised by a Legendre polynomial  No well established resonances are expected either in the pp or in the J/ψp and J/ψp channels. However, some resonances could potentially decay into pp [25], e.g. the f J (2220) [32] and the X(1835) [33,34]; thus they have been included in alternative models. The simplest model used to fit the data has no resonant contributions in the P + c , P − c and X decay chains, and is denoted as the baseline model. This model includes a nonresonant (NR) contribution in the X decay sequence with spin-parity quantum numbers equal to J P = 1 − , which has S-wave terms in both its production and decay. Indeed, due to the low Q-value of the decay, the S-wave contribution is expected to be favoured since higher values of orbital momentum are suppressed. Models including NR contributions with different quantum numbers (i.e. J P = 0 ± , 1 + ) are excluded because their −2 log L values are significantly worse than that of the J P = 1 − hypothesis.
Due to the limited sample size, the baseline model is described by two independent LS couplings for both B 0 s → J/ψX and X → pp decays, where L is the decay orbital angular momentum, and S is the sum of spins of the decay products. Fixing the two lowest orbital momentum couplings as the normalisation choice and three parameters, which are consistent with zero, reduces the number of free parameters to three.
The fit results of the baseline model are shown in Fig 2. The baseline model does not describe the data distribution, with a χ 2 goodness-of-fit test result of χ 2 /ndf = 64/38 corresponding to a p-value of 4 × 10 −5 . Therefore, two resonant contributions from P + c and P − c are added, with identical masses, widths and couplings. First, the P c (4312) state previously observed by the LHCb experiment in the Λ 0 b → J/ψpK − analysis [2] is included in the model with mass and width fixed at their known values. The broad P c structure with a mass around 4380 MeV, observed in 2015 [1], is not considered in this fit, since the helicity formalism used in Ref. [35] requires modifications in order to properly align the half-integer spin particles of different decay chains and, thus, those results need to be confirmed with an updated analysis of [36,37]. In this analysis no evidence for the P c (4312) state is found since the p-value, computed from the −2∆ log L of the alternative fit with respect to the default model, is measured to be 0.5. Exploiting the CL s method [38], an upper limit on the modulus of its coupling is set to 0.043 at 90% of confidence level, which corresponds to a fit fraction of 2.86%. A model with a new P ± c state given a free mass and width is chosen as the default model. Different spin-parity hypotheses for the P c states are investigated, i.e. J P = 1/2 ± and J P = 3/2 ± . Due to a limited sample size, only the lowest values of L are considered and the same coupling is assumed for all J P hypotheses, resulting in two free parameters: the modulus A(P c ) and the phase φ(P c ) of the coupling. The seven fit parameters ω contain the baseline model parameters, see Eq. 2, the coupling [A(P c ), φ(P c )], the mass and width of the P c state.
The fit result for the J P = 1/2 + hypothesis of the P + c state is shown in Fig. 2. The χ 2 /ndf is 36.7/36.8, where the number of degrees of freedom ndf is determined from fits to the χ 2 distribution extracted from pseudoexperiments. The statistical significance is estimated from pseudoexperiments generated with the baseline model and fitted with the default model, using amplitude parameters determined by the fit to data. The mass and width of the P c states are not defined in the baseline model, thus multiple fits to the same pseudodata are performed to account for the look-elsewhere effect, scanning the initial mass value in intervals of size 50 MeV. The test statistic t is built as the maximum of the −2 log L difference between the baseline and the default model [39] among all the fits obtained by scanning the initial mass values. The p-value is computed using a frequentist method as the fraction of pseudoexperiments with t larger than the t data value from the fits to data. The p-value ranges between 0.02% and 0.2% for different J P hypotheses, the lowest being associated to 1/2 + and the highest to 3/2 + , as reported in the Supplemental material. These p-values correspond to a signal significance in the range of 3.1 to 3.7σ, providing evidence for a new pentaquark-like state. Using the CL s method [38], none of the J P hypotheses considered can be excluded at 95% confidence level.
The hypothesis of a glueball state with mass equal to 2230 MeV and width of around 20 MeV [10] is also tested, by adding to the default model a resonance in the X decay chain with fixed mass and width. No evidence of f J (2220) is observed, as the fit with this contribution gives a p-value, computed from the −2∆ log L with respect to the default model, of 0.75 and an associated complex coupling of [−0.04 ± 0.09, −0.06 ± 0.16].
Systematic uncertainties are evaluated for the mass, width, coupling, and fit fractions of the sum of the P ± c contributions. For each source of uncertainty, pseudoexperiments are generated according to the alternative model with the same sample size as in data. The fit to such pseudoexperiments is performed using the default model. The systematic uncertainties, listed in Table 1, are assigned as the mean of the residual distributions between the fitted and the default parameter results. The main contributions are due to different NR models for the X decay chain, alternative J P hypotheses for the P c state and possible mismodelling of the efficiency distribution. The systematic uncertainty associated to the NR model is obtained including, in addition to the NR term with J P = 1 − and lowest values of L allowed, a P -wave resonant contribution with J P = 0 − , modelled with a Breit-Wigner lineshape in order to account for possible resonances, such as the X(1835) [33,34], decaying to a pp final state. Since none of the J P hypotheses investigated for the P ± c state can be excluded, an additional systematic uncertainty is assigned as the difference between the least and the most significant hypotheses. Finally, the uncertainty associated with the efficiency parameterisation is evaluated by summing two contributions. The first is obtained by replacing the default efficiency map with one determined from simulation of different data-taking conditions, and the second by using a parameterisation given by the product of one-dimensional functions of the considered fit variables. Other systematic uncertainties include alternative parameterisation of the background shape and the uncertainty in the background normalisation, which is varied within its statistical uncertainty. The background is parameterised using data in a sideband region around the B 0 s invariant-mass peak with m(J/ψpp) ∈ [5300, 5350] MeV and m(J/ψpp) ∈ [5420, 5460] MeV, to account for variations of the background as a function of the invariant mass. The default value of the hadron radius size for the Blatt-Weisskopf coefficients [40], equal to 3 GeV −1 , is replaced by two alternate values, 1.5 GeV −1 and 5 GeV −1 . Fit biases in the parameters estimation are extracted from the residual distribution of the generated and fitted parameters of pseudoexperiments based on the default model. Systematic uncertainties from orbital momentum for the NR, P c contributions, and invariant-mass resolution are found to be negligible. More details about systematic uncertainties can be found in the Supplemental material. The final significance including systematic uncertainties is equal to 3.1σ, which is the minimal value among the different sources of systematic uncertainty, as reported in Table 1.
The mass and width of this new pentaquark-like state are measured to be where the first uncertainty is statistical and the second systematic. The analysis of flavour untagged B 0 s decays is not sensitive to the P + c and P − c contributions separately, therefore −0.11 and phase φ(P c ) consistent with zero, corresponding to a fit fraction of (22.0 +8.5 −4.0 ± 8.6)% for the P c states. Due to the limited sample size, it is not possible to distinguish among different J P quantum numbers. A state compatible with this P c state is predicted in Ref. [41] with J P = 1/2 + .
In conclusion, an amplitude analysis of B 0 s → J/ψpp decays is presented, using data collected with the LHCb detector between 2011 and 2018, and corresponding to an integrated luminosity of 9 fb −1 . No evidence is seen for either a P c state at a mass of 4312 MeV [2] or the glueball state f J (2220) predicted in Ref. [10]. Unlike in other B decays [42][43][44][45], no threshold enhancement is observed in the pp invariant-mass spectrum, which is well modelled by a nonresonant contribution. Evidence for a Breit-Wigner shaped resonance in the J/ψp and J/ψp invariant masses is obtained with a statistical significance in the range of 3.1 to 3.7σ, depending on the assigned J P hypothesis.

Supplemental material A. Matrix element model
The amplitude model is constructed using the helicity formalism [29] where a two body decay A → (1)(2) contributes to the amplitude with a term where λ is the particle helicity, defined as the projection of the spin J A onto the momentum direction, and H A→(1)(2) are complex helicity couplings which describe the decay dynamics and take into account the Jacob-Wick phase factor of particle (2) [30] 2 and are defined as H . The Wigner D-matrix rotates the initial coordinate system of particle A, with the z axis aligned along the helicity axis of A, to the coordinate system with the z axis aligned along the particle 1 helicity axis. The angles φ 1 and θ 1 (known as "helicity angle" of A) represent the azimuthal and polar angles of particle 1 in the rest frame of A. The third angle is set to zero by convention. The last term of Eq. 4, R(m 12 ), is the lineshape dependence that contains either the Blatt-Weisskopf coefficients and threshold factors or, if A is a resonant contribution, a Breit-Wigner lineshape. Details are given below.
The helicity couplings are expressed in the LS basis, where L is the orbital and S is the spin angular momentum, using Clebsch-Gordan coefficients, B L,S , as For strong decays, possible L values are constrained by the conservation of parity: The four-body phase space of the B 0 s → J/ψpp decay is described by the pp invariant mass, m(pp), the two helicity angles θ p , θ µ of the p and µ in their parent reference frame and the azimuthal angle ϕ between the dihadron and dilepton decay planes.
In order to align different decay sequences, the cyclic ordering of the final particles is adopted to define the helicity angles, as suggested in Ref. [31], and the Jacob-Wick convention for particle (2). The ordering is important in order to guarantee a proper spin matching of the final particles and to ensure that all resonances share the same y axis instead of the opposite axis. For the decay under study, the following ordering is considered: where the angles are defined with respect to particle (1) for the P + c chain, to particle (2) for the X chain and to particle (3) for the P − c chain. In particular, in the X chain the p direction defines the X helicity angle, while in the P − c →pJ/ψ chain thep direction defines the P − c helicity angle. The Jacob-Wick phase-factor is also needed in order to properly align the final spin. Indeed, the rotation in the X → pp decay aligns the spin axis along the p momentum, while the rotation in the P + c → J/ψp frame aligns the spin axis in the direction opposite to the p momentum. Therefore, an additional 2 Particle (2) is the particle with opposite momentum with respect to particle (1) in the particle-A rest frame: p rotation to align the z axis between the P + c and X chains generates the particle (2) phase factor equal to (−1) Jp−λp in the amplitude of the P + c chain, where J p and λ p are the spin and the helicity of the proton in the P + c rest frame, respectively. Denoting the J/ψ as ψ, the matrix element for the B 0 s → XJ/ψ chain is where the first Wigner-D matrix of the B 0 s decay is omitted because the B 0 s has spin zero, implying that, for conservation of total angular momentum, the helicity of the ψ is equal to the helicity of X (|λ X − λ ψ | ≤ 0). The y axis of the ψ rest frame for all three decay sequences is chosen to be parallel to the y axis in the B 0 s rest frame to ensure a correct alignment. The angles φ are the polar and azimuthal angles of the µ − momentum in the ψ helicity frame. In the B 0 s rest frame, the proton momentum projected on the x axis is positively defined in order to satisfy φ p = 0. Since the ψ decay occurs through an electromagnetic interaction, the difference of the muon helicities, ∆λ, restricts to the values: ∆λ = λ µ − − λ µ + = ±1. As the value of ∆λ = 0 is suppressed by m µ /m ψ for the electromagnetic transition, it is omitted in the summation. The helicity coupling of the ψ decay is, therefore, ignored because it can be absorbed in the B 0 s couplings. The factor R(m 2 pp ) is the lineshape term of the pp invariant mass, as described below.
For J P (X) = 1 − , which is the best choice to fit the data, there are 3 (2) independent couplings H B 0 s λ ψ ,λ X ( H X λp,λ p ) to fit, already reduced by the parity conservation in the strong decay of the X resonance. Those couplings are then expressed in the LS basis; the lowest L and S of the B 0 s decay, B B 0 s L min S min , is always fixed to (1, 0) for each contribution. Therefore, the NR contribution contains two complex couplings for the B 0 s → J/ψX decay, associated to relative angular momenta L B 0 s = 0, 1, 2 and one for the decay X → pp, with L X = 0, 2, due to parity conservation. Here, L B 0 s and L X refer to the relative angular momentum between J/ψX and pp final states, respectively. In addition, due to the limited sample size, the couplings relative to L B 0 s = 1 for the production and to the phase φ of L X = 2 are consistent with zero and, therefore, are fixed to zero in the default model, reducing the number of free parameters to three. The number of LS couplings and free parameters are summarised in Table 2.
Similarly, the matrix elements for the P + c and P − c decay chains are given by where, as above, for angular momentum conservation, the following relations between helicities hold: λ P + c = λ p and λ P − c = λ p for the P + c and P − c chains, respectively. The Table 2: LS couplings and free parameters for the NR 1 −− model and the models with the P c states for different J P hypotheses. Here, (L, S) and (l, s) refer to the relative angular momentum and spin of P + c p (P − c p) and J/ψp (J/ψp) final states, respectively, for the P + c (P − c ) decay chain. The free parameters are denoted as (A, φ), where A is the modulus and φ the phase of the coupling.

Model
LS couplings x-axes are chosen in order to have φ are the polar and azimuthal angles of the µ − momentum in the ψ helicity frame in the two decay chains. Due to the limited sample size, from the fit to data neither two unique couplings for the P c states nor a relative phase between them can be extracted. A relation between the helicity couplings of P + c and P − c is thus imposed where the sign in the second line comes from the permutation of the final helicities, while the one of the first line depends on the weak dynamics and is chosen as such. This choice is arbitrary and does not impact the results because, due to limited sample size, there is no sensitivity to the relative interference between P ± c states. It is made such that it simplifies the model by constraining the two P ± c couplings to be equal and obtains the same interference pattern for all J P hypotheses, as discussed below in Eq. 23. The helicity couplings expressed in the LS basis are reduced to the lowest L values allowed. There is one independent coupling to fit, which is reduced to one free parameter for J P = 1/2 + and 3/2 − because the phase φ of the coupling is consistent with zero due to the limited sample size, as summarised in Table 2.
A single resonant contribution in the P c decay chains is parametrised by where is the Breit-Wigner function which includes the mass-dependent width Here, p is the momentum of the P c resonance in the B 0 s rest frame, q is the momentum of one of the decay products in the P c rest frame, while the momenta p 0 and q 0 denote their values at the resonance peak (m = M 0 ), and d is the hadron radius size fixed to 3 GeV −1 in the default fit and varied to 1.5 GeV −1 and 5 GeV −1 in the systematic uncertainty evaluation. The orbital momentum of the B 0 s decay is denoted as L B 0 s and that of the P c decay as L P cj . In order to properly account for the suppression due to higher values of L, the Blatt-Weisskopf coefficients are used together with the orbital barrier factor p L B L (p, p 0 , d). For NR contributions in the X chain, BW(m) is set to 1 and M 0(N R) to the midrange mass.
To sum the amplitudes from different decay chains coherently, the final state helicities of p, p and µ in the P c chains must be rotated in order to match the helicities of the X chain where θ P ± c p ( θ P ± c p ) are the polar angles in the p (p) rest frame between the boost directions from the P ± c to the X rest frames and α µ is the azimuthal angle in the ψ rest frame to correct the muon helicity states in the two chains. Regarding the definition of the polar angles, for the P + c chain, θ P + For the P − c chain, the θ P − c p and θ P − c p angles are obtained by substituting p ↔ p. Since the y axis is outgoing from the plane, the rotation to align the P + c decay chain to the X decay chain is counterclockwise between the J/ψ and thep momenta. While in thep rest frame, the rotation from the direction of P + c to that of p is clockwise, hence of an angle −θ P + c p . For the P − c chain, the rotation is always clockwise in the proton rest frame and counterclockwise in thep rest frame. For the J/ψ decay, since the muons come from the J/ψ for both decay chains, the polar angle is 0, implying the following: d 1/2 λ Pc µ ,λµ = δ λ Pc µ ,λµ . However, there is an azimuthal angle α µ because of the offset in the x axis. Since the boost to the µ rest frame is the same for both decay chains, i.e. always from the J/ψ rest frame, α µ can be determined in the J/ψ rest frame as where the index 3 refers to the rest frame after rotations;ẑ ψ 3 =p ψ µ and thex 3 axis can be derived asx as well asx The term aligning the muon helicity states between the two reference frames is given by The transformation for the µ − is equal to that of the µ + where the azimuthal angle takes a negative sign: α µ − = −α µ + . Considering the transformations of both muons, the final rotation is the multiplication of the two exponentials This implementation is found to be equivalent to the one proposed in Ref. [31], where the spin of the J/ψ in the different chains is aligned with a polar rotation before the J/ψ → µ + µ − decay.
Since the flavour of the B 0 s candidate is not tagged, the overall amplitude is the average of the B 0 s and B 0 s amplitudes, where the B 0 s amplitude is equal to that of B 0 s due to the absence of direct CP violation and is obtained inverting particles with antiparticles, i.e. p ↔ p, µ − ↔ µ + and inverting all azimuthal angles, φ ↔ −φ. To impose CP symmetry, the flavour eigenstates are projected onto the CP basis of B L and B H . Enforcing CP symmetry conservation, the B 0 s and B 0 s and their decay products can be written in terms of CP eigenstates. If CP is applied to the B 0 s final states, then CP |XJ/ψ = +|XJ/ψ , is obtained. Indeed, the contribution in the X chain is already a CP eigenstate, and in particular for J P = 1 − a CP -even one, while the two single P c are not. However, a combination of the P + c and P − c contributions can be projected onto a CP eigenstate orthogonal to the X contribution (CP -odd) by choosing the specific phase convention: e iφ = −1(−1) J Pc −3/2 , hence resulting in which spin dependent factor follows from Eq. 9. Therefore, the interference between the X and the combination of the P c contributions cancels out upon integration over the Dalitz plane, as they are orthogonal CP eigenstates.

B. Event-by-event efficiency parameterisation
Event-by-event acceptance corrections are applied to the data using an efficiency parameterisation based on the decay kinematics. The 4-body phase space of the topology B 0 s → J/ψ(→ µ + µ − )pp is fully described by four independent kinematic variables, one mass and three angles: m pp , θ p , θ µ , ϕ, where the angles are defined as, • θ µ and θ p : the helicity angles defined in the dimuon and dihadron rest frames, respectively; • ϕ: the azimuthal angle between the two decay planes of the dilepton and dihadron systems.
Since the final state is as self-conjugate, the p and the µ − are chosen to define the angles, for both B 0 s and B 0 s . For the signal mode, the overall efficiency, including trigger, detector acceptance and selection procedure, is obtained from simulation as a function of the four kinematic variables, ω ≡ {m pp , cos θ µ , cos θ p , ϕ }. Here, m pp and ϕ are transformed such that all four variables in ω lie in the range (−1, 1]. The efficiency is parameterised as the product of Legendre polynomials ε( ϕ) = i,j,k,l c i,j,k,l P (cos θ µ , i)P (cos θ p , j)P (ϕ , k)P (m pp , l), where P (x, n) is a Legendre polynomial of order n in x ∈ (−1, 1]. Employing the order of the polynomials as {3, 7, 7, 5} for {m pp , cos θ µ , cos θ p , ϕ }, respectively, was found to give a good parameterisation. The coefficients, c i,j,k,l , are determined from the simulation using a moments technique employing the orthogonality of Legendre polynomials where w n is the per-event weight taking into account both the generator level phase-space element, dΦ, and the kinematic event weights. Simulation samples are employed where B 0 s → J/ψpp events are generated uniformly in phase space. In order to render the simulation flat also in m(pp), the inverted phase-space factor, 1/dΦ, is considered. The factors of (2a + 1)/2 arise from the orthogonality of the Legendre polynomials, The sum in Eq. 24 is over the reconstructed events in the simulation sample after all selection criteria. The factor C ensures appropriate normalisation and it is computed such that where N rec is the total number of reconstructed signal events. Up to statistical fluctuations, the parameterisation follows the simulated data in all the distributions.

C. Significance studies for different J P hypotheses of the P c state
The significance studies for different J P hypotheses are reported in Table 3, together with the fit results. The significance is computed using a frequentist approach, by counting the number of pseudoexperiments above the −2∆ log L observed in data. The p-value and corresponding two-sided significance are reported in Tab. 3. Table 3: Results for all spin-parity hypotheses. Values of −2∆ log L (-2DLL), p-value and two-sided significance extracted from pseudoexperiments, together with mass (in MeV), width (in MeV), fit fraction and complex coupling (A,φ) of the P c states are reported.

D. Systematic uncertainties studies
In this section, more details about systematic uncertainties are presented. The systematic uncertainties are studied using pseudoexperiments generated according to the alternative model with the same statistics as in data. Parameters are then fitted with the default model. For each parameter, the uncertainty is taken as the mean value of the residual distribution between the fitted ensemble and the generated one. The main systematic uncertainties are described in the document. To estimate the uncertainty associated to the NR model parameterisation for the X decay chain, alternative NR models are used, also with different quantum numbers (J P = 0 ± , 1 + ). Combinations of NR and a resonant term are used to account for possible resonances decaying to a pp final state. The model, which yields comparable results in terms of − log L, comprises a NR term with J P = 1 + and a Breit-Wigner resonance with J P = 0 − , and is used to assign the systematic uncertainty. The impact of this model on the significance of the P c state is also studied and is found to be negligible. The other alternative models have been excluded based on the -2∆ log L, with respect to the baseline model, which is greater than 40. In addition, amplitude models including all values of orbital momentum for both the X and the P c decay chains are considered and their effect is found to be negligible. For the systematic uncertainty associated with the efficiency parameterisation, two different contributions are summed in quadrature. First, the default efficiency map extracted from all simulation events is replaced by using simulations from the Run 1 and the Run 2 data-taking periods, separately. Second, a simplified parameterisation with onedimensional Legendre polynomials is considered, neglecting the effect of correlations among the phase-space variables. Other sources of uncertainty regard different J P assignment for the P c state, the background parameterisation, the hadron radius size and the fit bias. The last uncertainty is estimated from pseudoexperiments by generating and fitting the parameters of the default model. Possible biases on the fit parameters are accounted for as the mean of the residual distributions between the generated and the fitted values.