Observation of $\psi(3686) \to n\bar{n}$ and improved measurement of $\psi(3686) \to p \bar{p}$

We observe the decay $\psi(3686) \to n \bar{n}$ for the first time and measure $\psi(3686) \to p \bar{p}$ with improved accuracy by using $1.07\times 10^8$ $\psi(3686)$ events collected with the BESIII detector. The measured branching fractions are $\mathcal{B}(\psi(3686) \to n \bar{n}) = (3.06 \pm 0.06 \pm 0.14)\times 10^{-4}$ and $\mathcal{B}(\psi(3686) \to p \bar{p}) = (3.05 \pm 0.02 \pm 0.12) \times 10^{-4}$. Here, the first uncertainties are statistical and the second ones systematic. With the hypothesis that the polar angular distributions of the neutron and proton in the center-of-mass system obey $1+\alpha \cos^2\theta$, we determine the $\alpha$ parameters to be $\alpha_{n\bar{n}} = 0.68 \pm 0.12 \pm 0.11$ and $\alpha_{p\bar{p}} = 1.03 \pm 0.06 \pm 0.03$ for $\psi(3686)\to n\bar{n}$ and $\psi(3686)\to p\bar{p}$, respectively.


I. INTRODUCTION
As a theory of the strong interaction, QCD has been well tested in the high energy region. However, in the lower energy region, nonperturbative effects are dominant, and theoretical calculations are very complicated. The charmonium resonance ψ(3686) has a mass in the transition region between the perturbative and nonperturbative regimes. Therefore, studying ψ(3686) hadronic and electromagnetic decays will provide knowledge of its structure and may shed light on perturbative and nonperturbative strong interactions in this energy region [1]. Nearly four decades after the decay ψ(3686) → pp was measured [2], we are able to measure ψ(3686) → nn for the first time using the large ψ(3686) samples collected at BESIII [3]. A measurement of ψ(3686) → nn, along with ψ(3686) → pp, allows the testing of symmetries, such as flavor SU(3) [4].
Due to the Okubo-Zweig-Iizuka mechanism, the decays of J/ψ and ψ(3686) to hadrons are mediated via three gluons or a single photon at the leading order. Perturbative QCD predicts the "12% rule," [36,37]. This rule is expected to hold for both inclusive and exclusive processes but was first observed to be violated in the decay of ψ into ρπ by MARKII [38], called the "ρπ puzzle." Reviews of the relevant theoretical and experimental results [39][40][41] conclude that the current theoretical explanations are unsatisfactory. Further precise measurements of J/ψ and ψ(3686) decay to NN may provide additional knowledge to help understand the ρπ puzzle.
In this paper, we report the first measurement of ψ(3686) → nn and an improved measurement of ψ(3686) → pp. First, we introduce the BESIII detector and the data samples used in our analysis. Then, we describe the analysis and results of the measurements of ψ(3686) → nn and ψ(3686) → pp. Finally, we compare the branching fractions and α values with previous experimental results and different theoretical models.

II. BESIII DETECTOR, DATA SAMPLES AND SIMULATION
BESIII is a general purpose spectrometer with 93% of 4π solid angle geometrical acceptance [42]. A small cell, helium-based multilayer drift chamber (MDC) provides momentum measurements of charged particles with a resolution of 0.5% at 1 GeV/c in a 1.0 T magnetic field and energy loss (dE/dx) measurements with a resolution better than 6% for electrons from Bhabha scattering. A CsI(Tl) electromagnetic calorimeter (EMC) measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end caps). A time-of-flight system (TOF), composed of plastic scintillators, with a time resolution of 80 ps (110 ps) in the barrel (end caps) is used for particle identification (PID). A superconductive magnet provides a 1.0 T magnetic field in the central region. A resistive plate chamber-based muon counter located in the iron flux return of the magnet provides 2 cm position resolution and is used to identify muons with momentum greater than 0.5 GeV/c. More details of the detector can be found in Ref. [42]. This analysis is based on a ψ(3686) data sample corresponding to 1.07×10 8 events [3] collected with the BESIII detector operating at the BEPCII collider. An offresonance data sample with an integrated luminosity of 44 pb −1 [3], taken at the c.m. energy of 3.65 GeV, is used to determine the non-ψ(3686) backgrounds, i.e. those from nonresonant processes, cosmic rays, and beamrelated background.
A Monte Carlo (MC) simulated "inclusive" ψ(3686) sample of 1.07 × 10 8 events is used to study the background. The ψ(3686) resonance is produced by the event generator kkmc [43], while the decays are generated by evtgen [44,45] for the known decays with the branching fractions from the particle data group [26], or by lundcharm [46] for the remaining unknown decays. Signal MC samples for ψ(3686) → NN are generated with an angular distribution of 1 + α cos 2 θ, using the α values obtained from this analysis. The interaction of particles in the detectors is simulated by a geant4-based [47] MC simulation software boost [48], in which detector resolutions and time-dependent beam-related backgrounds are incorporated.

III. MEASUREMENT OF ψ(3686) → nn
The final state of the decay ψ(3686) → nn consists of a neutron and an antineutron, which are back to back in the c.m. system and interact with the EMC. The antineutron is expected to have higher interaction probability and larger deposited energy in the EMC. To suppress background efficiently and keep high efficiency for the signal, a root-based [49] multivariate analysis (MVA) [50] is used.

A. Event selection
A signal candidate is required to have no charged tracks reconstructed in the MDC. Events are selected using information from the EMC. Showers must have deposited energy of E > 25 MeV in the barrel (| cos θ| < 0.8) or E > 50 MeV in the end caps (0.86 < | cos θ| < 0.92). The "first shower" is the most energetic shower in the EMC, and the first shower group (SG) includes all showers within a 0.9 rad cone around the first shower. The direction of a SG is taken as the energy-weighted average of the directions of all showers within the SG. The SG's energy, number of crystal hits and moments are the sums over all included showers for the relevant variables. The "second shower" is the next most energetic shower excluding the showers in the first SG, and the second SG is defined based on the second shower analogous to how the first SG is defined. The "remaining showers" are the rest of the showers which are not included in the two leading SGs.
We require | cos θ| < 0.8 for both SGs, and the energies of the first SG and second SG to be larger than 600 MeV and 60 MeV, respectively. The larger energy requirement applied to the first SG is to select the antineutron, which is expected to have larger energy deposits in the EMC than the neutron due to the annihilation of the antineutron in the detector. There is a total of 6 × 2 + 2 = 14 variables, which are listed in Table II, including the energies, number of hits, second moments, lateral moments, numbers of showers, largest opening angles of any two showers within an SG, and number and summed energy of the remaining showers.
We implement the MVA by applying the boosted decision tree (BDT) [51]. Here, 50 × 10 3 signal and 100 × 10 3 background events are used as training samples. The signal events are from signal MC simulation, and the background events are a weighted mix of selected events from the off-resonance data at √ s = 3.65 GeV, inclusive MC simulation, and exclusive MC simulation samples of the processes ψ(3686) → γχ cJ , χ cJ → nn, (J = 0, 1, 2), which are not included in the inclusive MC samples. The scale factors are 3.7 for the off-resonance data, determined based on luminosity and cross sections [3], and 1.0 for the inclusive MC sample. We also select independent test samples with the same components and number of events as the training samples. The "MVA" selection Here, r0 = 5 cm is the average distance between crystal centers in the EMC, ri is the radial distance of crystal i from the cluster center, and Ei is the crystal energy in decreasing order.

B. Background determination
The signal will accumulate in the large θ open region since the final states are back to back. The possible peaking background of e + e − → γγ is studied with a MC sample of 10 6 events. After the final selection, and scaled to the luminosity of real data, only 27±10 events are expected from this background source, which can be neglected. This is also verified by studying the off-resonance data. The remaining backgrounds are described by three components, which are the same as those used in the BDT training. None of them produces a peak in the θ open distribution.

C. Efficiency correction
The neutron and antineutron efficiencies are corrected as a function of cos θ in the e + e − c.m. system to account for the difference between data and MC simulation. Control samples of ψ(3686) → pnπ − + c.c., selected using charged tracks only, are used to study this difference.
The efficiency of the BDT selector for the antineutron is defined as ǫ = N BDT /N tot , where N tot is the total number of antineutron events obtained by a fit to the pπ recoil mass distribution, and N BDT is the number of antineutrons selected with the BDT method. The same shower variables as used in the nominal event selection are used in the BDT method to select the antineutron candidate. The efficiency for the neutron is determined analogously. The ratios of the efficiencies of MC simulation and data as a function of cos θ are assigned as the correction factors for the MC efficiency of the neutron and antineutron, and are used to correct the event selection efficiencies. The ratios and corrected efficiencies are shown in Fig. 1 for the neutron and antineutron separately. The corrected efficiencies are fitted by fourthorder polynomial functions with χ 2 /ndf = 0.87 and 1.13 for the neutron and antineutron, respectively.

D. Branching fraction and angular distribution
We perform a fit to the θ open distribution of data to obtain the numbers of signal candidates and background events. The histogram from signal MC simulation is used to construct the signal probability density function (PDF). Corresponding histograms from the three background components, as described in Sec. III B, are used to construct the background PDFs. The numbers of events from each source are free parameters in the fit. Figure 2 shows the fit to the θ open distribution. The fit yields N sig = 6056 ± 117 nn events with χ 2 /ndf = 3.24. Using a corrected efficiency ǫ = 18.5%, the branching fraction of ψ(3686) → nn is determined to be (3.06 ± 0.06) × 10 −4 via B = N sig /(N ψ(3686) ǫ), where N ψ(3686) is the total number of ψ(3686) and the uncertainty is statistical only.
We fit the cos θ n and cos θn distributions separately with fixed fractions of each component to determine the α values. For these fits, an additional selection criterion θ open > 3.01 is used to further suppress the continuum background, and the fractions of each components within the region 3.01 < θ open < 3.20 are obtained from the θ open fit results. For the cos θ n and cos θn distributions, the background PDFs are constructed with the same method as used in the fits to θ open , while the signal PDF is constructed by the formula (1 + α cos 2 θ)ǫ(θ). Here, ǫ(θ) is the corrected polar angle-dependent efficiency parameterized in a fourth-order polynomial, as described in Sec. III C. Figure 3 shows the fits to the cos θ n and cos θn distributions. An average α nn = 0.68 ± 0.12 for the angular distribution is obtained, while the separate fit results are 0.76 ± 0.12 (χ 2 /ndf = 0.81) and 0.60 ± 0.12 (χ 2 /ndf = 2.01) for the cos θ n and cos θn distributions, respectively. The uncertainties here are statistical only. Since the neutron and antineutron are back to back in the c.m. system and the two angular distributions are fully correlated, the average does not increase the statistics, and the uncertainty is not changed.

Resolution of θopen
To determine the difference in the θ open resolution between data and MC, we fit the θ open distribution of data with the signal PDF convolved with a Gaussian function of which the parameters are left free in the fit. The resultant mean and width of the Gaussian function are 0.005 and 0.002 rad, respectively. With these modified PDFs, the resultant changes are 0.3% for the branching fraction and 0.0% for the α value, which are taken as the systematic uncertainties from the resolution of θ open . We do not consider the resolution effect for the cos θ distributions because of their smoother shapes.

Backgrounds
The uncertainties associated with the background amplitudes are estimated by fitting the θ open distribution with fixed contributions for the continuum and inclusive MC background. The differences between the new results and the nominal ones, 0.8% and 8.1% for the branching fraction and the α value, respectively, are taken as the systematic uncertainties related with the background amplitudes.
To estimate the effect on the α distribution from the continuum background shape, we redo the fit to the cos θ distributions with the shape of the continuum background obtained without the additional requirement θ open > 3.01, assuming that there is no correlation between θ open and cos θ. The difference of α to the nominal result is 4.4%.
All in all, we determine the uncertainty from backgrounds to be 0.8% for the branching fraction and 9.2%, the quadratic sum of 8.1% and 4.4%, for α.

Neutral reconstruction efficiencies
The reconstruction efficiency is corrected in bins of cos θ, and the uncertainty of the correction is taken to be the statistical uncertainty, which is about 2% per cos θ bin. To obtain its effect on our results, we allow the efficiency to fluctuate about the corrected efficiency according to the statistical uncertainty, and redo the fits with the modified efficiencies. We also use the histograms of the corrected MC efficiencies directly. The largest change of the signal yield is 0.2% with the average efficiency changing by 2% (1% each fromn and n), and the largest change in α is 12.8%. We take these differences from the standard results as the systematic uncertainties of the neutral efficiency correction.  The fit result is shown as the solid blue curve. The signal shape is from MC simulation and is presented as the dashed black histogram. The background is described by three components: continuum background in dotted red, inclusive MC sample in dash-dotted green, and the tiny contribution from ψ ′ → γχcJ , χcJ → nn (not included in the inclusive MC sample) in long-dashed cyan. All yields are free parameters in the fit.

Remaining showers
We have checked and found that the number and energy of remaining showers are independent of the angle, as we expected. Then only the branching fraction measurement will be affected by the unperfect MC simulation. Based on the distributions of the number and energy of remaining showers from the data, we weight them in the signal MC considering their correlation. The difference is found to be 0.4% by comparing the efficiencies obtained with and without weighting, and is quoted as the corresponding uncertainty.

Analysis method
We perform input/output checks by generating different signal MC samples with different α values, from zero to unity; mixing these signal MC samples with backgrounds; and scaling these samples to the number of events according to data. Compared to the input values, the output signal yield is very close to the input, and its corresponding systematic uncertainty can be neglected. For the measurement of α, the average difference, 2%, is taken as the systematic uncertainty. The fit result is shown as the solid blue curve. The signal shape is parametrized by (1 + α cos 2 θ)ǫ(θ), shown as the dashed black curve. The background is described by three components: continuum background in dotted red, inclusive MC sample in dash-dotted green, and a small contribution from ψ ′ → γχcJ , χcJ → nn in long-dashed cyan.

Binning
In the nominal analysis, the θ open , cos θ n , and cos θn distributions are divided into 60, 20 and 20 bins, respectively. To estimate the uncertainty associated with binning, we redivide the distributions of θ open , cos θ n , and cos θn into [55, 65], [18,22], [18,22] bins, respectively and perform 11 × 5 × 5 = 275 fits of θ open , cos θ n , and cos θn with all possible combinations of binnings to determine the signal yields and α values. The differences between the average results and the nominal values, 0.1% for the branching fraction and 4.5% for the α value, are taken as the systematic uncertainties.

Physics model
The signal efficiency in the branching fraction measurement depends on the value of α. Varying α by its standard deviation, the relative change on the detection efficiency, 1.1%, is taken as the systematic uncertainty due to the physics model.

Trigger efficiency
The neutral events used for this analysis are selected during data taking by two trigger conditions: 1) the number of clusters in the EMC is required to be greater than one, and 2) the total energy deposited in the EMC is greater than 0.5 GeV [52]. The efficiency of the former condition is very high [52], and we conservatively take 2% as its systematic uncertainty. Requiring the EMC total energy to be larger than 0.9 GeV, the trigger efficiency of the second condition is 98.8% [52], with an uncertainty of 1.2%. Comparing the nominal results to the results with the higher total energy requirement, the difference is 0.2%. Combining the two gives 1.4%, which is taken as the systematic uncertainty of the second trigger condition. Since these two trigger conditions may be highly correlated, we take a conservative 3.4% as the total systematic uncertainty of the trigger.

Summary of systematic uncertainties
The systematic uncertainties in the measurements of ψ(3686) → nn are summarized in Table III. Assuming these systematic uncertainties are independent of each other, the total uncertainty is obtained by adding the individual uncertainties quadratically.

A. Event selection
The final state of ψ(3686) → pp consists of a proton and an antiproton, which are back to back and with a fixed momentum in the c.m. system. A candidate charged track, reconstructed in the MDC, is required to satisfy V r < 1.0 cm and |V z | < 10.0 cm, where V r and V z are the distances of closest approach of the reconstructed track to the interaction point, projected in a plane transverse to the beam and along the beam direction, respectively. Two charged track candidates with net charge zero are required. We also require the momentum of each track to satisfy 1.546 < p < 1.628 GeV/c in the c.m. system, which is within three times the resolution of the expected momentum, and the polar angle to satisfy | cos θ| < 0.8. Using the information from the barrel TOF, likelihoods L i for different particle hypotheses are calculated, and the likelihood of both the proton and antiproton must satisfy L p > 0.001 and L p > L K , where L p is the PID likelihood for the proton or antiproton hypothesis, and L K is the likelihood for the kaon hypothesis. Further, we require the opening angle of the two tracks to satisfy θ open > 3.1 rad in the ψ(3686) c.m. system. There are 18, 984 candidate events satisfying the selection criteria, which are used for the further study.

B. Background estimation
In the analysis, backgrounds from the continuum process e + e − → pp and ψ(3686) decay into non-pp final states are explored with different approaches. The former background is studied with the off-resonance data at √ s = 3.65 GeV. With the same selection criteria, there are (22 ± 5) events that survive, and the expected background in the ψ(3686) data is (22 ± 5) × 3.7 = 81 ± 18 events, where 3.7 is the scale factor which is the same as in the nn study. By imposing the same selection criteria on the ψ(3686) inclusive MC sample, no non-pp final state events survive, and the non-pp final state background from ψ(3686) decays is negligible. We also check the latter background with the two-dimensional sidebands of the proton versus antiproton momenta, which is shown in Fig. 4. There are a few events in the sideband regions, marked as A and B in Fig. 4, but MC studies indicate that the events are dominantly initial state or final state radiation events of ψ(3686) → pp. The ratios of events in each sideband region to that in signal region are consistent between data and signal MC simulation.

C. Efficiency correction
In the ψ(3686) → pp analysis, we correct the MC efficiency as a function of cos θ of the proton and antiproton, where the corrected factors include both for tracking and PID efficiencies. The efficiency differences between data and MC simulation, which are obtained by studying the same control sample of ψ(3686) → pp, are taken as the correction factors. To determine the efficiency for the proton, we count the number of ψ(3686) → pp events by requiring an antiproton only, and then check if the other track is reconstructed successfully in the recoiling side and passes the PID selection criterion. The efficiency is defined as n 2 /(n 1 + n 2 ), where n 1 and n 2 are the yields of events with only one reconstructed track identified as an antiproton and with two reconstructed tracks identified as proton and antiproton, respectively. The yields n 1 and n 2 are obtained from fits to the antiproton momentum distributions. In the fit, the signal shape is described by the momentum distribution of the antiproton with the standard selection criteria for ψ(3686) → pp, and the background is described by a first-order polynomial function since it is found to be flat from a study of the inclusive MC sample. Cosmic rays and beamrelated backgrounds are subtracted using V z -sidebands, in which |V z | ≤ 5 cm is defined as the signal region and (−10 < V z < −5) and (5 < V z < 10) are defined as sideband regions. A similar analysis is performed for the antiproton detection efficiency. The ratio of efficiencies between MC simulation and data are displayed individually in Fig. 5 for the proton and antiproton. We obtain the corrected MC efficiency to select ψ(3686) → pp candidates, also shown in Fig. 5. The corrected MC efficiencies are fitted with fourth-order polynomial functions with χ 2 /ndf = 2.56 and 2.57 for the proton and antiproton, respectively.

D. Branching fraction and angular distribution
After subtracting the continuum background, the branching fraction is determined to be B(ψ(3686) → pp) = (3.05 ± 0.02) × 10 −4 via B = N sig /(N ψ(3686) ǫ) with the corrected efficiency of ǫ = 58.1% determined with the angular distribution corresponding to the value of α obtained in this analysis. The cos θ distributions of the proton and anti-proton for the selected candidates are shown in Fig. 6. The distributions are fitted with the functional form N sig (1 + α cos 2 θ)ǫ(θ) + N bg f bg , where N bg and f bg , the yield and the shape of the continuum background, are fixed in the fit according to the off-resonance data at    for the data and MC simulation, respectively, where the efficiencies are estimated by integrating the Gaussian function within the specific signal regions. Thus, the systematic uncertainty is taken to be 0.4% for the two charged tracks.
The total systematic uncertainty associated with the momentum resolution for the branching fraction is 0.5%, and that for the α value measurement is found to be negligible.

Background
The dominant background is from the continuum process, which is estimated with the off-resonance data sample at √ s = 3.65 GeV. The corresponding uncertainty of 18 events, which is 0.1% of all signal events, is taken as the uncertainty in the branching fraction measurement associated with the background. The uncertainty on the α value associated with background is studied by leaving the background yield free in the fit and found to be negligible.

Tracking and PID efficiencies
In the nominal analysis, the tracking and PID efficiencies for the proton and anti-proton are corrected to improve the accuracy of the measurement. Thus, only the uncertainties associated with the statistics of correction factors and the method to exact correction factors are considered.
We repeat the analysis 1000 times by randomly fluctuating the correction factors for the proton and antiproton detection efficiency with Gaussian functions independently in the different cos θ bins, where the width of the Gaussian function is the statistical uncertainty of the correction factors. The standard deviations of the results are < 0.1% for the branching fraction and 0.2% for α, which are taken as the systematic uncertainties associated with the statistical uncertainties.
In the nominal analysis, the corrected efficiency is parametrized with a fourth-order polynomial function. Alternative parametrizations with a polynomial function symmetric in cos θ and directly using the histogram for the corrected efficiency are performed. The maximum changes of the branching fraction and α value, 3.3% and 2.1%, respectively, are taken as the systematic uncertainties.
To be conservative, the linear sums of the two uncertainties, 3.3% and 2.3%, are taken as the systematic uncertainties for the branching fraction and α measurements associated with the tracking and PID efficiency, respectively.

Method
From input/output checks, the average relative differences between measured and true values are 1.1% for the branching fraction and 2.0% for α, which are taken as the systematic uncertainties.

Binning
In the nominal analysis, the cos θ range of the proton and anti-proton of (−0.8, 0.8) is divided into 20 bins to determine the corrected tracking and PID efficiency. Alternative analyses with 10 or 40 bins are also performed, and the largest differences with respect to the nominal results are taken as the systematic uncertainties associated with binning. The effect is negligible for the branching fraction measurement and 1.0% for the α measurement. · · · · · · Number of ψ(3686) 0.7 · · · Total 4.0 3.2

Physics model
In the branching fraction measurement, the detection efficiency depends on the value of α. Alternative detection efficiencies varying α from 0.96 to 1.10, corresponding to one standard deviation, are used. The largest change of the efficiency with respect to the nominal value, 1.8%, is taken as the systematic uncertainty.

Trigger efficiency
Events with two high momentum charged tracks in the barrel region of the MDC have trigger efficiencies of 100.0% and 99.94% for Bhabha and dimuon events [52], respectively, and the systematic uncertainty from the trigger is negligible.

Summary of systematic uncertainties
The systematic uncertainties of ψ(3686) → pp from the different sources are summarized in Table IV. Assuming the systematic uncertainties are independent, the total uncertainty is the sum on the individual values added in quadrature.

V. SUMMARY AND DISCUSSION
In this paper, we measure the branching fractions of ψ(3686) → nn and pp, and the α values of the polar angle distribution, which are described by 1+α cos 2 θ. The final results are B(ψ(3686) → nn) = (3.06±0.06±0.14)×10 −4 and α nn = 0.68 ± 0.12 ± 0.11, and B(ψ(3686) → pp) = (3.05 ± 0.02 ± 0.12) × 10 −4 and α pp = 1.03 ± 0.06 ± 0.03, where the former process is measured for the first time and the latter one has improved precision compared to previous measurements, as summarized in Table I. The measured α pp is close to 1.0, which is larger than previous measurements, but both B(ψ(3686) → pp) and α pp are consistent with previous results within the uncertainties.
To check for an odd cos θ contribution from the 2γ exchange process [53], we fit the angular distributions as before but with the function 1 + β cos θ + α cos 2 θ. The results are β nn = 0.04 ± 0.05 and β pp = 0.01 ± 0.02. The possible contributions from odd cos θ terms in this analysis are consistent with zero.
With the assumption the decay process is via a single photon exchange, the α value must satisfy |α| ≤ 1 [54]. Then, the formula 1+sin φ cos 2 θ is applied to fit to the pp data again, and we obtain the result φ pp = 1.57 ± 0.28 ± 0.25, where the statistical uncertainty is obtained from fit directly and the systematical uncertainty is propagated from the 3.2% of the α pp value.
In the decay of J/ψ → nn and pp [26], both the branching fractions and α values are very close between the two decay modes, which is expected if the strong interaction is dominant in J/ψ → NN decay and the relative phase of between the strong and electromagnetic amplitudes is close to 90 • [13]. In contrast, in ψ(3686) decays, the branching fractions are quite close between the two decay modes, but the α values are not, which may imply a more complex mechanism in the decay of ψ(3686) → NN . It makes a similar and straightforward extraction of the phase angle impossible in the decay of ψ(3686) → NN , and further studies are deserved.