Production of $b\bar{b}$ at forward rapidity in $p$+$p$ collisions at $\sqrt{s}=510$ GeV

The cross section of bottom quark-antiquark ($b\bar{b}$) production in $p$+$p$ collisions at $\sqrt{s}=510$ GeV is measured with the PHENIX detector at the Relativistic Heavy Ion Collider. The results are based on the yield of high mass, like-sign muon pairs measured within the PHENIX muon arm acceptance ($1.2<|y|<2.2$). The $b\bar{b}$ signal is extracted from like-sign dimuons by utilizing the unique properties of neutral $B$ meson oscillation. We report a differential cross section of $d\sigma_{b\bar{b}\rightarrow \mu^\pm\mu^\pm}/dy = 0.16 \pm 0.01~(\mbox{stat}) \pm 0.02~(\mbox{syst}) \pm 0.02~(\mbox{global})$ nb for like-sign muons in the rapidity and $p_T$ ranges $1.2<|y|<2.2$ and $p_T>1$ GeV/$c$, and dimuon mass of 5--10 GeV/$c^2$. The extrapolated total cross section at this energy for $b\bar{b}$ production is $13.1 \pm 0.6~(\mbox{stat}) \pm 1.5~(\mbox{syst}) \pm 2.7~(\mbox{global})~\mu$b. The total cross section is compared to a perturbative quantum chromodynamics calculation and is consistent within uncertainties. The azimuthal opening angle between muon pairs from $b\bar{b}$ decays and their $p_T$ distributions are compared to distributions generated using {\sc ps pythia 6}, which includes next-to-leading order processes. The azimuthal correlations and pair $p_T$ distribution are not very well described by {\sc pythia} calculations, but are still consistent within uncertainties. Flavor creation and flavor excitation subprocesses are favored over gluon splitting.

The cross section of bottom quark-antiquark (bb) production in p+p collisions at √ s = 510 GeV is measured with the PHENIX detector at the Relativistic Heavy Ion Collider. The results are based on the yield of high mass, like-sign muon pairs measured within the PHENIX muon arm acceptance (1.2 < |y| < 2.2). The bb signal is extracted from like-sign dimuons by utilizing the unique properties of neutral B meson oscillation. We report a differential cross section of dσ bb→µ ± µ ± /dy = 0.16 ± 0.01 (stat) ± 0.02 (syst) ± 0.02 (global) nb for like-sign muons in the rapidity and pT ranges 1.2 < |y| < 2.2 and pT > 1 GeV/c, and dimuon mass of 5-10 GeV/c 2 . The extrapolated total cross section at this energy for bb production is 13.1 ± 0.6 (stat) ± 1.5 (syst) ± 2.7 (global) µb. The total cross section is compared to a perturbative quantum chromodynamics calculation and is consistent within uncertainties. The azimuthal opening angle between muon pairs from bb decays and their pT distributions are compared to distributions generated using ps pythia6, which includes nextto-leading order processes. The azimuthal correlations and pair pT distribution are not very well described by pythia calculations, but are still consistent within uncertainties. Flavor creation and flavor excitation subprocesses are favored over gluon splitting.

I. INTRODUCTION
The bottom-quark production in hadron-hadron collisions is an important test of perturbative quantum chromodynamics (pQCD) calculations. Because of its large mass, m b Λ QCD , the b-quark production cross section can be accurately calculated by including next-toleading order (NLO) processes, especially at high center of mass energies [1]. The measurement of the bb production cross section over a wide range of colliding energies in hadron-hadron collisions provides a stringent test of pQCD theory calculations and a baseline measurement for studying modifications of heavy quark production in heavy ion collisions.
Cross section measurements for bottom production in hadron-hadron collision experiments have been made from lower energy fixed-target experiments [2][3][4] ( √ s < 45 GeV) up to collider energies ( √ s > 100 GeV). It was found that pQCD predictions match experimental results well at energies greater than √ s = 1 TeV [5][6][7][8][9][10][11][12], but less so at lower energies. Results at the wide range of collision energies of the Relativistic Heavy Ion Collider explore an important gap between the low-energy fixed-target and TeV-energy regimes.
Without displaced vertex b-tagging capability at PHENIX, b-quark production has been studied using unlike-sign dileptons from heavy quark decays [13]. The PHENIX and STAR collaborations have previously measured the bottom cross section in p+p collisions at √ s = 200 GeV using electron-hadron correlations [14,15] and using dilepton invariant mass and momentum distributions [16][17][18].
Like-sign dimuons have previously been used to investigate the phenomenon of neutral B meson oscillations at e + e − colliders by the CLEO Collaboration [19], the AR-GUS Collaboration [20], the ALEPH Collaboration [21] * Deceased † PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov and the UA1 Collaboration [22] in p +p collisions. In this measurement, we use the yield of like-sign dimuons along with the properties of neutral B meson oscillation to determine the bb cross section.
In the Standard Model, neutral B meson oscillation is a result of higher order weak interactions that transform a neutral B meson into its antiparticle: B 0 →B 0 because the flavor eigenstates differ from the physical mass eigenstates of the meson-antimeson system [23,24]. In the absence of oscillation, as shown in Fig. 1(a), primary-primary decays, where the lepton's direct parent is the B meson, can only produce unlike-sign lepton pairs. For example b →B(B − ,B 0 ,B 0 s , ..) → l − and b → B(B + , B 0 , B 0 s , ..) → l + while like-sign lepton pairs can result from a mixture of primary and secondary decays (decay chain).
However, if oscillation occurs, as is the case for neutral B mesons (B 0 d and B 0 s ), theB 0 meson can spontaneously change into a B 0 meson as shown in Fig. 1(b). Unless otherwise noted, we denote B(B) as a generic admixture of bottom (antibottom) hadrons with production ratios, from weak decays (i.e. Z → bb) of: B + (B − ) = 40.4±0.9%, B 0 (B 0 ) = 40.4±0.9%, B 0 s (B 0 s ) = 10.3±0.9%, and b(b)-baryon = 8.9 ± 1.5%. The B c production ratio is negligible (0.2%) and less than the uncertainties associated with bottom hadrons listed above [25]. The timeintegrated probability for a neutral B meson to oscillate before it decays is defined as where ∆m is the mass difference between heavy and light mass eigenstates and Γ is the decay rate of the weak eigenstates. These values are found to be χ d ≈ 0.1874 ± 0.0018 and χ s ≈ 0.499311 ± 0.000007 for the B 0 d and B 0 s mesons, respectively [25]. This process can result in a like-sign dilepton event from a primary-primary decay (see Fig. 1(b)). Given the large branching ratio of the B → µ decay channel (≈ 10.99%) [25], the like-sign dilepton from a primary-primary decay provides a unique opportunity for extracting the bb cross section.
In this paper, we present measurements of bb production cross section through the like-sign dimuon decays and the azimuthal opening angle between the muon pair and their p T distributions in p+p collisions at √ s = 510 GeV at forward (1.2 < y < 2.2) and backward (−2.2 < y < −1.2) rapidities. The azimuthal opening angle and pair p T distributions are compared to distributions generated using pythia6 with parton-shower (ps) model [26]. The model approximates the correction to all higher orders (almost next-to-leading-log) for bb production, which includes flavor creation, flavor excitation, and gluon splitting. The extrapolated total cross section, using ps pythia6 [26] and pythia8 [27], and MC@NLO [28] calculations, is also presented and compared to pQCD calculation.
The paper is organized as follows: The PHENIX apparatus is described in Sec. II. The data samples used for this analysis and the analysis procedure are presented in Sec. III. The results are presented and discussed in Sec. IV. The summary and conclusions are presented in Sec. V

II. EXPERIMENTAL SETUP
A complete description of the PHENIX detector can be found in Ref. [29]. We briefly describe here only the detector subsystems used in these measurements. The relevant systems, which are shown in Fig. 2, include the PHENIX muon spectrometers covering forward and backward rapidities and the full azimuth. Each muon spectrometer comprises a hadronic absorber, a magnet, a muon tracker (MuTr), and a muon identifier (MuID). The absorbers comprise layers of copper, iron, and stainless steel and have about 7.2 interactions lengths. Following the absorber in each muon arm is the MuTr, which comprises three stations of cathode strip chambers in a radial magnetic field with an integrated bending power of 0.8 T·m. The MuID comprises five alternating steel absorbers and Iarocci tubes. The composite momentum resolution, δp/p, of particles in the analyzed momentum range is about 5%, independent of momentum and dominated by multiple scattering. Muon candidates are identified by reconstructed tracks in the muon spectrometers.
Another detector system relevant to this analysis is the beam-beam counter (BBC), consisting of two arrays of 64Čerenkov counters, located on both sides of the interaction point and covering the pseudorapidity range 3.1 < |η| < 3.9. The BBC system was used to measure the p+p collision vertex position along the beam axis (z vtx ), with 2 cm resolution, and initial collision time. It was also used to measure the beam luminosity and form a minimum bias trigger (MB). The MB trigger requires at least one hit in each BBC on the sides of the interaction point.

A. Data set and quality cuts
The data set for this analysis is collected by PHENIX during the 2013 p+p run at √ s = 510 GeV. Events, in coincidence with the MB trigger, containing a muon pair within the acceptance of the spectrometer are selected by the level-1 dimuon trigger (MuIDLL1-2D) requiring that at least two tracks penetrate through the MuID to its last two layers. After applying a vertex cut of |z vtx | < 30 cm and extensive quality assurance checks, the data remaining correspond to 3.02 × 10 12 MB events or to an integrated luminosity of 94.4 pb −1 .
A set of cuts was used to select good muon candidates and improve the signal-to-background ratio. Hits in the MuTr are used to make MuTr tracks and hits in the MuID are used to make MuID roads. The MuTr track is required to have more than 9 hits out of the maximum possible of 16 while the MuID road is required to have more than 6 hits out of the maximum possible of 10. Additional χ 2 cut is applied on MuTr track that is calculated from the difference between the measured hit positions of the track and the subsequent fit. MuTr tracks are then projected to the MuID at the first MuID gap and matched to MuID roads by applying cuts on maximum position and angle differences. Muon candidates are required to have a minimum p T greater than 1 GeV/c. This cut improves the sample quality by reducing background from pions and kaons. A minimum of 3.0 GeV/c is applied to single muon momentum along the beam axis, p z , which is reconstructed and energy-loss corrected at the collision vertex, corresponding to the momentum cut effectively imposed by the absorbers. Muon candidates are further restricted to the rapidity range of −2.2 < y < −1.2 for the south muon arm and 1.2 < y < 2.2 for the north muon arm. Additionally, a cut on the χ 2 of the fit of the two muon tracks to the common vertex of the two candidate tracks near the interaction point is applied.

B. Detector acceptance and reconstruction efficiency
The product of the acceptance and reconstruction efficiency (A ) is determined using Monte Carlo (MC) simulation. The A is defined by the number of dimuons reconstructed in the muon spectrometers with respect to the number of dimuons generated in the same kinematic region. The kinematic distributions of pythia 1 [30] generated p T , rapidity, and bb mass shape were used as input into a full PHENIX geant4 simulation [31].
The p T and rapidity distributions were tuned such that the reconstructed distributions match those of 2013 data. Variations within the uncertainties of data are taken as systematic uncertainty.
The detector response in the simulation is tuned to a set of characteristics (dead and hot channel maps, gains, noise, etc.) that describes the performance of each detector subsystem. The simulated vertex distribution is also tuned to match that of the 2013 data. The simulated events are further embedded with real data to account for the effects of detector noise and other background tracks, and then are reconstructed in the same manner as the real data. The A as a function of dimuon mass is 1 We used pythia6 (ver 6.421), with parton distribution functions given by CTEQ6LL. The following parameters were modified: A as a function of invariant mass for like-sign dimuons (a), as a function of dimuon azimuthal opening angle (b) and as a function of dimuon pT (c). These distributions are the weighted average of µ + µ + and µ − µ − distributions.
shown in Fig. 3 (a), as a function of dimuon opening angle is shown in Fig. 3 (b) and as a function of dimuon p T is shown in Fig. 3 (c). The relative difference in A between the two spectrometers is due to different detection efficiencies of the MuTr and MuID systems and different amount of absorber material.

C. Raw yield extraction
We measure like-sign dimuons in the same muon arm that have an invariant mass between 5 and 10 GeV/c 2 . In this mass range, the dimuon spectrum includes correlated and uncorrelated pairs. The correlated pairs are dominated by the semi-leptonic decay of open bottom pairs from the primary-secondary decay chain (see Fig. 1(a)) or from the primary-primary pairs from neutral B meson oscillation (see Fig. 1(b)). Dileptons from the Drell-Yan process and quarkonia decays can only yield unlike-sign pairs. D mesons can produce likesign pairs through their decay chain. For example, c → D + → µ + + anything and the other open charm decays asc → D − → K + + anything → µ + ν µ . However, in the mass range of interest the like-sign pairs from D mesons are negligible. The contribution from neutral D meson oscillation to the like-sign signal is expected to be very small because the oscillation probability is O(< 10 −2 ) [32]; therefore, it is not included.

Correlated background
Additional contribution to the correlated pairs could originate from correlated sources such as dijets or punch through hadrons. Hadrons (particularly π ± and K ± ) can punch through to the last gap of the MuID or decay to muons creating a background to the correlated like-sign signal. These contributions are estimated using MC simulation by determining the p T -dependent survival probability that a hadron will traverse the muon arm detectors and applying it to pythia generated dihadron pairs to get the yield expected at the back of the muon arm detectors. π ± and K ± are generated with pythia 2 [18,30,33] and then run through the PHENIX detector simulation chain to determine a p T -dependent probability that the hadrons penetrate the last gap of the MuID.
To get a better estimate of the survival probability, the hadron simulation is run using two different hadron interaction packages for geant: fluka and geisha [34,35]. Figure 4 shows the simulated invariant mass spectra from irreducible background are fitted with an exponential function of the form exp(a + b × m + c × m 2 ) between 5 and 10 GeV/c 2 , where m is the invariant mass and a, b and c are fit parameters. The average of the indicated results from geisha and fluka is used to subtract the hadronic background from like-sign pairs while the difference is considered as a systematic uncertainty.
The invariant mass distribution for like-sign pairs is then constructed from pythia generated dihadron pairs within the same event and from mixed events, with each entry weighted by the survival probability. Event-mixing procedure is discussed in the next section. Just as with data, the correlated like-sign signal is obtained by subtracting the mixed event spectrum from the like-sign spectrum, providing the correlated like-sign signal due to dijets or punch through hadrons. The sum of π and K correlated like-sign signals is weighted based on their p T -dependent cross sections [36,37].
Fake like-sign pairs due to charge misidentification and like-sign pairs from Drell-Yan process or quarkonia decays and muon-decayed or punch through hadrons were also studied and found to be negligible.

Uncorrelated background
The uncorrelated pair contribution is estimated using event mixing technique [38], where like-sign pairs are constructed by pairing muons in the current event with those of the same sign and same arm in previous events of z-vertex position within 2 cm. The mixed event pairs (N BG ++ and N BG −− ) form the uncorrelated background spectrum which is normalized to the the foreground (N F G ++ and N F G −− ) using a normalization factor The normalization factor requires that the integrated counts from event mixing equal those from the like-sign in the low mass region where the correlated pairs are expected to be negligible [38]. The normalized like-sign pairs from event mix-ing are given as: However, the specific range where the signal of interest is negligible is not well known, and the average of normalization factors over five mass ranges (0.6-2.6 GeV/c 2 , 1.0-2.0 GeV/c 2 , 1.6-3.2 GeV/c 2 , 2.6-3.8 GeV/c 2 , and 0.6-4.2 GeV/c 2 ) is used. The correlated like-sign signal (N cor ±± ) is then isolated by subtracting the mixed-event spectrum (N BG ±± ) from the "foreground" like-sign pairs (N F G ±± ) according to the following, To further improve the normalization process, the bb invariant mass distribution shape from ps pythia6 simulation is utilized. This is done by normalizing the integral of the ps pythia6 distribution to the result of Eq. (3), over the signal mass range 5-10 GeV/c 2 . The integral of the normalized bb mass distribution is then subtracted from the background distribution in Eq. (2) for each of the background ranges and the normalization factor is recalculated. The second step is then repeated until the value of the mixed-events normalization factor converges. Figure 5 shows the resulting distributions of N F G ±± , N BG ±± and N cor ±± as a function of the invariant mass of the pairs. These distributions are corrected with A . To extract the bb distribution as a function of the azimuthal opening angle between muon pairs (∆φ) and their p T , the normalization factors obtained previously are used to normalize ∆φ and p T mixed event distributions, which are then subtracted from ∆φ and p T foreground distributions, respectively. Table I summarizes the systematic uncertainties. Evaluated as standard deviations, they are divided into three categories based upon the effect each source has on the measured results:

D. Systematic uncertainties
Type-A: Point-to-point uncorrelated uncertainties that allow the data points to move independently with respect to one another and are added in quadrature with statistical uncertainties; however, no systematic uncertainties of this type are associated with this measurement.
Type-B: Point-to-point correlated uncertainties which allow the data points to move coherently within the quoted range to some degree. These systematic uncertainties include a 4% uncertainty from MuID tube efficiency and an 8.2% (2.8%) from MuTr overall efficiency for the north (south) arm. The systematic uncertainty associated with A includes the uncertainties on the input p T and rapidity distributions which are extracted by varying these distributions over the range of the statistical uncertainty of the data, yielding 4.4% (5.0%) for the north (south) arm. To be consistent with the real data analysis, a trigger emulator was used to match the MuIDLL1-2D trigger for the data. The efficiency of the trigger emulator was studied by comparing the dimuon mass spectrum requiring the dimuon passes the trigger emulator to the dimuon mass spectrum requiring the dimuon passes the MuIDLL1-2D trigger, which resulted in a 1.5% (2%) uncertainty for the north (south) arm. Additional 11.2% (8.8%) systematic effect for the north (south) arm was also considered to account for the azimuthal angle distribution difference between data and simulation.
The source of systematic uncertainty in signal extraction is the normalization of mixed events which could come from the choice of the different normalization ranges in the mixed events or bb shape from pythia used to guide the signal extraction. A 1.9% uncertainty on the mixed events normalization was observed from using each of the five normalization windows by itself as well as the different combinations of these normalization windows. pythia bb shape is the sum of three subprocesses: flavor creation, flavor excitation and gluon splitting. A maximum variation of 3.1% on the extracted signal was observed from choosing each of the subprocesses by itself as the source of bb shape. Added in quadrature, they result in a 3.6% uncertainty on signal extraction.
The systematic uncertainty associated with correlated backgrounds could come from the input p T distribution, differences between geisha and fluka, and differences between geant3 and geant4. pythia p T distributions of π ± and K ± were compared separately to fits of UA1 data [36,37] and an overall difference of 18% was observed. Differences of up to 30% and 20% between fluka and geisha, see Fig. 4, were observed in the north and south arms, respectively. Additional 15% was considered to account for the difference between geant3 and geant4. Added in quadrature, all three sources give an overall effect on the hadronic background of 39% (31%) for the north (south) arm for the mass and ∆φ distributions. For p T distribution, a p T -dependent correction was used for the effect on the input p T spectra and the other two sources gave an overall effect on the hadronic background of 34% (25%) for the north (south) arm. To extract the systematic uncertainty associated with the cross section (or invariant yields) for all distributions (mass, ∆φ and p T ), the hadronic background was varied between the limits listed above which resulted in an overall systematic of 5.1% (4.5%) for north (south) arm.
The Type-B systematic uncertainties are added in quadrature and amount to 16.0% (12.8%) for the north (south) arm. They are shown as shaded bands on the associated data points.
Type-C: An overall (global) normalization uncertainty of 10% was assigned for the BBC cross section and efficiency uncertainties [39] which allows the data points to move together by a common multiplicative factor.

IV. RESULTS AND DISCUSSION
A. Differential cross section The differential yield and cross section of B meson pairs decaying into like-sign dimuons as a function of where N µµ /A (y, m) is the yield of like-sign dimuons from B meson decay normalized by A (y, m) in y and m bin with ∆y and ∆m widths, respectively. The yields of the north and south arms are calculated independently and are consistent within statistical uncertainties; therefore, the weighted average [40] is used in the differential yield calculation. MB BBC = 0.53 ± 0.02 is the fraction of inelastic p+p collisions recorded by the BBC. σ pp BBC = 32.5 ± 3.2 mb is the cross section as seen for the BBC in p+p collisions at √ s = 510 GeV, which is determined from the van der Meer scan technique [41]. BBC = 0.91 ± 0.04 is the efficiency of the MB trigger for events containing a hard scattering [42]. N MB is the number of MB events.
The differential cross section of like-sign dimuons from B meson decay is shown in Fig. 6. The gray shaded bands represent the weighted average of the quadratic sum of type-B systematic uncertainties of the north and south arms, ≈10.1%. The average is weighted based on the statistical uncertainties of each arm. In addition to type-B systematic uncertainties, we have a 10% global systematic uncertainty for BBC cross section and efficiencies [39].
To obtain the differential cross section of all B meson pairs that decay into dimuons, regardless of the muon pair charge, the differential cross section of likesign dimuons from B meson decay is scaled by the ratio of the total number of all B meson pairs that decay into dimuons, regardless of their sign, to those of like-sign. For clarification purposes, the process is divided into two separate steps defined by two variables α(m) and β, both of which depend on the signal from like-sign dimuons due to oscillation.
The ratio of like-sign dimuons at mass m and from primary-primary decays due to B 0 oscillation to like-sign muon pairs resulting from primary-primary or a mixture of primary-secondary decays is defined as: which is calculated in the mass range 5 < m < 10 GeV/c 2 at 1.2 < |y| < 2.2 and p T > 1 GeV/c and extrapolates the correlated like-sign signal to an open bottom signal from oscillation, N osc ±± . The α(m) is obtained using open bottom events from three model calculations: mc@nlo (ver 4.10), ps pythia6 (ver 6.421) and pythia8 (ver 8.205), as shown in Fig. 7. The red line is a secondorder polynomial fit with χ 2 /ndf of 3.8/4. The shaded boxes represent the uncertainty based on the three model calculations.
β is the ratio of primary-primary like-sign dimuons due to B 0 oscillation to all B meson pairs that decay into primary-primary dimuons with all possible muon charge pairs (++, −− and +−). β converts the number of muon pairs from oscillation into all B meson pairs and is defined as: The value of β is 0.22 ± 0.01 which is the calculated RMS value from the three model simulations described above. The error of β is the standard deviation of the three model calculations which represents the modeldependent uncertainty.  Differential cross section of all dimuons from B meson decay. The error bars represent the statistical uncertainties, and the gray shaded band represents the quadratic sum of type-B systematic uncertainties.
The differential cross section of all B meson pairs that decay into a primary-primary dimuon, regardless of the muon pair charge, is then calculated as follows: ditional type-B systematic uncertainties associated with this measurement due to α(m) and β and amount to 1.9% and 4.5%, respectively, are included. This brings the type-B systematic uncertainties on d 2 σ bb→BB→µµ /dydm to 11.2%.

B. Total cross section
To extrapolate from the bb differential cross section in the muon decay channel within the acceptance of muon arms to a total bb cross section, the differential cross section is scaled by the ratio of B pairs that decay to dimuons within the measured region to those over the entire kinematic range. This method is similar to that found in Ref. [44].
The total cross section, σ bb , is extrapolated and corrected for the semi-leptonic branching ratio in the following manner: where BR B→µ is the branching ratio of B to muon through the primary decay channel (=10.99%), and scale, defined as: which is a factor used to convert from the visible kinematic region to full phase space. The scale factor is determined from pythia and mc@nlo simulations. It is taken as the average value, 1.96 × 10 −3 , of ps pythia6 (CTEQ6LL), ps pythia6 (CTEQ5M1), pythia8 (CTEQ6LL) and mc@nlo (CTEQ5M), as listed in Table II. The difference in the scale factor due to the different models and parton distribution functions is considered to be a global type-C uncertainty, which amounts  to 18.1%. This results in a total cross section of 13.1 ± 0.6 (stat)±1.5 (type-B syst)±2.7 (global syst) µb. Type-B systematic uncertainties are from the differential cross section while global uncertainties are the quadrature sum of type-C from the differential cross section and uncertainties arising from the extrapolation.
The σ bb measured at √ s = 510 GeV is shown in Fig. 9(a) and compared to measurements from other experiments [2-4, 9, 17, 45]. The solid line is the cross section from NLO pQCD calculations [43] and the dashed lines are error bands, and they are obtained by varying the renormalization scale, factorization scale and bottom quark mass. At √ s = 510 GeV, the NLO pQCD calculation predicts σ bb = 11.5 +6.5 −3.9 µb, which is consistent with the extrapolated total cross section using the current dimuon analysis within uncertainties. Figure 9(b) shows the ratio of data to theory.

C. Azimuthal correlations and pair pT
The like-sign µµ pair yield from bb decays is shown in Fig. 10 and Fig. 11 as a function of ∆φ and pair p T , respectively. The spectra are compared to model calculations based on ps pythia6 that are normalized by fitting the subprocesses sum to the data [18]. The generated pairs are filtered with the same kinematic cuts that are applied in the data analysis.
The azimuthal opening angle distribution from ps pythia6 shows a similar pattern to that of the data, an increase until ≈2.6 rad and then drop, and it is consis-tent with the data with χ 2 /ndf ≈27.3/28. The data show steeper p T dependence than that of ps pythia6 but they are still consistent when considering the large statistical and systematic uncertainties. These results show similar behavior to that observed at 200 GeV [18] where the data favors a dominant mix of flavor creation and flavor excitation subprocesses over gluon splitting.

V. SUMMARY AND CONCLUSION
In summary, we presented first measurements of the differential cross section for dimuons from bottom quarkantiquark production in p+p collisions at √ s = 510 GeV, which we found to be: dσ bb→µ ± µ ± /dy = 0.16 ± 0.01 (stat) ± 0.02 (syst) ± 0.02 (global) nb. The analysis technique is based on the yield of high-mass correlated like-sign dimuons from parity-violating decays of B meson pairs at forward and backward rapidities. Using a model dependent extrapolation, the measured differential cross section is converted to a total cross section of 13.1±0.6 (stat)±1.5 (syst)±2.7 (global) µb. This extrapolated total cross section is consistent with NLO pQCD calculations within uncertainties.
The azimuthal opening angle between the muons from bb decays and the pair p T distributions are compared to distributions generated using ps pythia6 [26], which includes NLO processes. While the data tend to have a wider azimuthal distribution than ps pythia6 calculations and present a steeper p T distribution, both are still consistent within uncertainties with ps pythia6, where flavor creation and flavor excitation subprocesses are dominant. This is similar to what was observed at 200 GeV [18].