1 Introduction

Measurements of Semi-Inclusive Deep-Inelastic Scattering (SIDIS)

$$\begin{aligned} {\mu }+{N}\rightarrow \mu '+ nh + X, \quad n=1,2, \ldots \end{aligned}$$
(1)

of high-energy polarised muons \(\mu \) off nucleons N in the initial state and scattered muons \(\mu '\), n measured hadrons h and unobserved particles X in the final state are sensitive to the spin-dependent Parton Distribution Functions (PDFs) of nucleons. The SIDIS cross section depends, in particular, on the azimuthal angle of each produced and measured hadron (see e.g. Refs. [1,2,3,4,5]), which leads to azimuthal asymmetries related to convolutions of the nucleon Transverse-Momentum-Dependent (TMD) PDFs and parton-to-hadron Fragmentation Functions (FFs). These asymmetries can appear in SIDIS off unpolarised, longitudinally or transversely polarised nucleons.

The TMD PDFs were studied in a number of experiments. The short overview of earlier results obtained by the HERMES, CLAS and COMPASS collaborations on azimuthal asymmetries in SIDIS production of charged hadrons was given in Refs. [6, 7]. The COMPASS collaboration has published results on asymmetries off unpolarised \(^6\)LiD (referred to as “deuteron”) target [10], transversely polarised deuterons [11, 12] and transversely polarised NH\(_3\) (referred to as “proton”) target [13,14,15,16,17,18]. The common analysis of transversely polarised deuteron and proton data is included in Refs. [13,14,15,16,17,18] also. The updated overview of the TMD PDFs including the COMPASS results can be found in Ref. [19]. The COMPASS results on azimuthal asymmetries off longitudinally polarised deuterons based on the data collected in 2002, 2003 and 2004 were published in Ref. [6] for so called “integrated” asymmetries and asymmetries as functions of kinematic variables extracted in the restricted kinematic region. Similar data have been collected in 2006 but not published yet. The results on the integrated asymmetries for 2006 and for the combined 2002–2006 data which are presented in this Paper are extracted using hadrons from the whole available kinematic region, at variance with Ref. [6, 7], while asymmetries as functions of the kinematic variables are extracted using the hadrons from a restricted kinematic region, similar to Ref. [6].

The Paper is organised as follows. The SIDIS kinematics, basic formulae and a brief theoretical overview are given in Sect. 2. The analysis of the 2006 data is described in Sect. 3. The results on the asymmetries of the combined 2002–2006 data are presented in Sect. 4. Systematic uncertainties are discussed in Sect. 5 and conclusions are given in Sect. 6.

2 Theoretical framework

The SIDIS kinematics is illustrated in Fig. 1.

Fig. 1
figure 1

The SIDIS kinematics shown for target deuteron polarisation \({\varvec{P}}_\Vert \) antiparallel to the beam direction

The 4-momenta of the incident and scattered muons are denoted by l and \(l'\), respectively. The 4-momentum of the virtual photon is given by \(q=l-l'\) with \(Q^2=-q^2\). The angle of the momentum vector \(\varvec{q}\) of the virtual photon with respect to the incident muon is denoted by \(\theta _\gamma \). The vectors \(\varvec{p}^h\) and \(\varvec{P}_\Vert \) denote the hadron momentum and the longitudinal target deuteron polarisation, respectively. Their transverse components \({\varvec{p}}_T^h\) and \({\varvec{P}_T}\) are defined with respect to the virtual-photon momentum. The longitudinal component \(|{\varvec{P}_L}|= |{\varvec{P}_\Vert }|\cos \theta _\gamma \) is approximately equal to \(|\varvec{P_\Vert }|\) due to the smallness of the angle \(\theta _\gamma \). The small transverse component is equal to \(\left| {\varvec{P}_T}\right| =|{\varvec{P}_{\Vert }}|\sin \theta _\gamma \) where \(\sin \theta _\gamma \approx 2({Mx}/{Q})\sqrt{1 - y}\). Here, M is the nucleon mass and \(y= (qp)/(pl)\) is the fractional energy of the virtual photon, where p is the 4-momentum of the target nucleon. The angle \(\phi \) denotes the azimuthal angle between the lepton scattering plane and the hadron production plane, while \(\phi _S\) denotes the angle of the deuteron polarisation vector with respect to the scattering plane: \(\phi _S=0^\circ \) or \(180^\circ \) for deuteron polarisation parallel or antiparallel to the beam direction, respectively. Furthermore, the Bjorken variable, \(x_{Bj}\equiv x =Q^2/(2pq)\), the fraction of the virtual-photon energy taken by a hadron, \(z=(pp^h)/(pq)\), the transverse momentum of a hadron, \(p^h_T\), and the invariant mass of the photon-nucleon system, \(W^2=(p+q)^2\), that, together with \(Q^2 > 1\, (\hbox {GeV}/c)^2\) and \(0<y< 1\), characterise SIDIS under study.

The general expression for the differential SIDIS cross section (see Refs. [1,2,3,4,5] and references therein) is a linear function of the incident muon polarisation \(P_\mu \) and of the longitudinal and transverse components \({\varvec{P}}_L\) and \({\varvec{P}}_T\) of the target deuteron polarisation \(\varvec{P}_\Vert \):

$$\begin{aligned} \mathrm {d}\sigma= & {} \mathrm {d}\sigma _{00}+{\varvec{P}_\mu } \mathrm {d}\sigma _{L0}+{\varvec{P}_L}\left( {\mathrm {d}\sigma _{0L} +{\varvec{P}_\mu } \mathrm {d}\sigma _{LL}}\right) \nonumber \\&+{\varvec{P}_T} \left( {\mathrm {d}\sigma _{0T}+{\varvec{P}_\mu } \mathrm {d}\sigma _{LT}}\right) . \end{aligned}$$
(2)

Here, the first (second) subscript of the partial cross sections refers to the beam (target) polarisation: \(0,\ L\) or T denote unpolarised, longitudinally or transversely polarised.

The azimuthal asymmetries of charged hadron production \(a_{h^\pm }(\phi )\) are defined as follows:

$$\begin{aligned} a_{h^\pm }(\phi )= \frac{\mathrm {d}\sigma ^{\leftarrow \Rightarrow }- \mathrm {d}\sigma ^{\leftarrow \Leftarrow }}{|P_L|(\mathrm {d}\sigma ^{\leftarrow \Rightarrow }+\mathrm {d}\sigma ^{\leftarrow \Leftarrow })}~, \end{aligned}$$
(3)

where all cross sections are functions of the angle \(\phi \). The Eq. (3) represents a definition of the experimentally measured asymmetries common for this Paper and for Refs. [6,7,8,9]. The first (second) superscript denotes the beam (target) spin orientation. The symbol \(\leftarrow \) denotes the incident muon spin orientation that, in the case of a positive charge of the incident muons, is mainly opposite to the beam direction. For the CERN muon beam, the average value of \(|{\varvec{P}_\mu }|\) is equal to 0.8. The beam polarisation does not enter in the definition of measured asymmetries. The symbols \(\Rightarrow \) and \(\Leftarrow \) denote the target deuteron spin orientations (polarisations) the first of which is parallel, considered further as positive \((+)\), and the second one is antiparallel \((-)\) to the beam direction (see Sect. 3.1).

Substituting the general expression for \(\mathrm {d}\sigma \) (Eq. 2) in the cross sections of the Eq. (3), one can obtain the expected contributions of the partial cross sections to the azimuthal asymmetries. As the result, when taking into account the signs of the target polarisations, one can see that only four partial cross sections contribute to the numerator of Eq. (3) and two to its denominator. In the numerator, we expect to have contributions from \(\mathrm {d}\sigma _{0L}\), \(P_\mu \mathrm {d}\sigma _{LL}\) and \(\tan \theta _\gamma (\mathrm {d}\sigma _{0T}+P_\mu \mathrm {d}\sigma _{LT})\), while in the denominator from \(\mathrm {d}\sigma _{00}\) and \(P_\mu \mathrm {d}\sigma _{L0}\). The explicit expression for these partial cross sections in terms of the PDFs and their dependences on the hadron azimuthal angle have been given in Refs. [6, 7] and briefly commented below.

Following phenomenological considerations based on the QCD parton model of the nucleon and SIDIS in one-photon exchange approximation, the squared modulus of the matrix element, defining the cross sections, is represented by a number of diagrams. As an exampleFootnote 1, the diagram accounting for the contribution to the SIDIS cross section of the chiral-odd transversity PDF \(h_{1}(x)\) convoluted with the chiral-odd Collins FF \(H_1^\bot (z)\) is shown in Fig. 2. The diagram contributes to the asymmetry Eq. (3) via the term \(\mathrm {d}\sigma _{0T}\). Other PDFs convoluted with corresponding FFs can also contribute to the cross section of spinless or unpolarised hadron production off longitudinally polarised deuterons and their expected azimuthal modulations. As motivated in Refs. [6, 7], out of predicted terms the \(\phi \)-independent term \(a^0_{h^\pm }\) and four modulation terms up to the order of M / Q are retained for the analysis:

$$\begin{aligned} a_{h^\pm }(\phi )= & {} a_{h^\pm }^{0} + a_{h^\pm }^{\sin \phi }\sin \phi + a_{h^\pm }^{\sin 2\phi }\sin 2\phi + a_{h^\pm }^{\sin 3\phi }\sin 3\phi \nonumber \\&+\, a_{h^\pm }^{\cos \phi }\cos \phi . \end{aligned}$$
(4)

The sign and the amplitude of each modulation is a subject of the \(a_{h^\pm }(\phi )\) data analysis (see Sect. 3.6). In Eq. (4), the term \(a_{h^\pm }^{0}\) is related to the well known helicity PDF \(g_{1L}\) contributing to asymmetries via \(d\sigma _{LL}\). The two terms with amplitudes \(a_{h^\pm }^{\sin 2\phi }\) and \(a_{h^\pm }^{\sin \phi }\) reported in Ref. [22] are related to the worm-gear-L PDF \(h_{1L}^\bot \), and to the PDFs \(h_L\) and \(f_L\), respectively, contributing to the asymmetries via \(d\sigma _{0L}\). The transversity PDF \(h_1\) and Sivers PDF \(f_{1T}^\bot \) can also contribute to the term with amplitude \(a_{h^\pm }^{\sin \phi }\) via \(d\sigma _{0T}\) with a small factor \(\tan \theta _\gamma \). Other terms in Eq. (4), not seen yet in experiments with longitudinally polarised deuterons, are the terms with amplitudes \(a_{h^\pm }^{\sin 3\phi }\) and \(a_{h^\pm }^{\cos \phi }\). They are related to the pretzelosity PDF \(h_{1T}^\bot \) and the worm-gear-T PDF \(g_{1T}\) contributing to asymmetries via \(d\sigma _{0T}\) and \(d\sigma _{LT}\), respectively, suppressed by the small factor \(\tan \theta _\gamma \). The COMPASS results [6] obtained from the 2002 – 2004 data suggested some indications for a possible x-dependence of terms with amplitudes \(a_{h^\pm }^{\sin 2\phi }\) and \(a_{h^\pm }^{\cos \phi }\). The contribution of the term with amplitude \(a_{h^\pm }^{\cos 2\phi }\) which could have appeared from \(\mathrm {d}\sigma _{00}\) in the denominator of Eq. (3) is disregarded. This amplitude is expected [23, 24] to be of the order of 0.1 and would enter Eq. (4) with the factor \(a^0_{h^\pm }\), that is of the order of \(10^{-3}\) for integrated asymmetries (see Table 2), or with \(a_{h^\pm }^0 (x)\le 0.05\) for asymmetries as functions of kinematic variables (see Fig. 6). This is beyond our experimental accuracy. The same comments apply to possible contributions of terms with amplitudes \(a_{h^\pm }^{\cos \phi }\) and \(a_{h^\pm }^{\sin \phi }\) which also could originate from \(\mathrm {d}\sigma _{00}\) and \(\mathrm {d}\sigma _{L0}\) of the denominator of Eq. (3), respectively. The negligible impact of the disregarded modulations on the amplitudes in Eq. (4) is confirmed by the 2006 data (see Sect. 3.6). All modulation amplitudes obtained in this Paper refer to the average value of the beam polarisation equal to −0.8.

Fig. 2
figure 2

The diagram describing the contribution to the SIDIS cross section of PDF \(h_1(x)\) convoluted with FF \(H_1^\bot (z)\)

The aim of this study is to continue searches for possible modulations in \(a_{h^\pm }(\phi )\) as manifestation of TMD PDFs describing the nucleons in the deuteron and to investigate the x, z and \(p_T^h\) dependences of the corresponding modulation amplitudes. For these purposes, we used first the 2006 deuteron data and then the combination of all 2002–2006 COMPASS deuteron data with longitudinal target polarisation.

3 The 2006 data analysis

3.1 Experimental set-up

The COMPASS set-up is a two-stage forward spectrometer with the world’s largest polarised target and various types of tracking and particle identification detectors (PID) in front and behind of two large-aperture magnets SM1 and SM2. These detectors provide data for reconstruction of corresponding tracks. The spectrometer was operated in the high energy (160 GeV) muon beam at CERN. Its initial configuration (see Ref. [20]) was used for data taking in 2002–2004. During the long accelerator shutdown in 2005, the set-up was modified (see Ref. [21]). The major modifications influencing the present analysis were as follows: (i) the replacement of the two 60 cm long target cells (denoted as U and D) by three cells \(U,\ M\) and D of lengths 30 cm, 60 cm and 30 cm, (ii) the replacement of the target solenoid magnet by the new one with a wider aperture and (iii) the installation of the electromagnetic calorimeter ECAL1 in front of the hadron calorimeter HCAL1. The ECAL1 is not used in the analysis because it was not fully operational yet in 2006 and partially acted as a hadron absorber. These modifications of the apparatus where aimed at further reduction of systematic uncertainties, enlargement of the spectrometer acceptance and improvement of \(e/\gamma \)  PID capabilities. These modifications have required reconsidering the Refs. [6, 7] methods of data stability tests and asymmetry calculations (see Sects. 3.3, 3.6.)

Table 1 Intervals of x, z, \(p^h_T\) and their weighted mean values for which asymmetries as functions of kinematic variables were calculated. The \(Q^2\)-intervals corresponding to the x-intervals are shown for reference

The data in 2006 were taken in two groups of periods. Each group is characterised by its initial set of polarisations in the target cells which are obtained by using different frequencies of the microwave field to polarise the target material (deuterons) in different cells at the certain direction of the target magnet solenoid field. The solenoid fields holds the polarisation. The field direction is denoted as \(f=+\), if it coincides with the beam direction, or \(f=-\), if opposite. The first group of the periods is denoted by G1 and other one by G2. Each period includes a certain number of intervals of continuous data taking (referred to as runs). The G1 data taking periods started with the initial setting of positive deuteron polarisation in target cells U and D and the negative one in cell M, both corresponding to \(f = +\). After taking some number of runs, the field was reversed to \(f=-\) causing the reversal of the target cell polarisations, so that the data were taken with opposite deuteron polarisations in the cells. The periodic reversal of polarisations continued up to the end of G1 periods. Within the periods, the cell polarisations, needed for asymmetry calculations (see Sect. 3.5), were measured for each run in order to make sure that they are stable at the level of about 55%. If polarisations dropped below this limit, they were restored by the microwave field before the beginning of the next period. For G2 periods, the procedure was analogous but the initial setting of polarisation in the cells was opposite to the one in G1 at the same field \(f=+\). The periodic reversal of the cell polarisations within each group of periods was used to estimate a possible time-dependent systematics of the data. The change of the initial setting of the cell polarisations was used to estimate a possible systematic change of the spectrometer acceptance due to superposition of the solenoid field and the field of SM1. If there is no such systematic change, the acceptance in G1 and G2 periods must be the same for stable performance of the spectrometer.

3.2 Selection of SIDIS events and hadrons

Let us call as “SIDIS event” an event determined by Eq. (1) and reconstructed with tracks using the data recorded by the tracking and PID detectors.

The overall statistics of 2006 is about \(44.6\times 10^6\) of preselected candidates for inclusive DIS and SIDIS event with \(Q^2>1~(\hbox {GeV}/c)^2\). The sample was obtained after rejection of runs that did not pass the data stability tests (see Sect. 3.3) and events that did not pass the reconstruction tests. The latter ones were rejected if the Z-coordinate (along the beam) of the interaction point (vertex) was determined with an uncertainty larger than \(3\sigma \) of average which varied within 1.5–2 cm for different target cells.

Fig. 3
figure 3

Kinematic distributions of selected SIDIS events vs. \(Q^2\) and y and of charged hadrons vs. z and \(p^h_T\) within the region shown in Table 1: 2006 (lower, red), 2002 – 2004 (middle, green) and 2002–2006 (upper, blue)

The selection of SIDIS events from the preselected sample was done as described in Refs. [6, 7]. For each SIDIS event, a reconstructed vertex with incident (\(\mu \)) and scattered (\(\mu ^\prime \)) muons and one or more additional tracks were required. Trajectories of the incident muons were required to traverse all target cells in order to have the same beam intensity for each of them. The track crossing more than 30 radiative lengths along the reconstructed trajectory was associated with \(\mu ^\prime \). The cuts were applied on the quality of the reconstructed tracks forming vertices, the effective lengths of the target cells (28 cm, 56 cm, 28 cm), the momentum of incident muons (\(140~\hbox {GeV}/c- 180~\hbox {GeV}/c\)), the fractional energy carried by all tracks from the event \((z < 1)\) and the fractional virtual-photon energy \((0.1< y < 0.9)\). About \(36.6\times 10^6\) SIDIS event candidates remained after cuts.

The distribution of track multiplicities per SIDIS candidates peaks at four. These tracks include scattered \(\mu ^\prime \) and hadron candidates. For a track to be identified as hadron, it was required that: its transverse momentum was larger than \(0.05~\hbox {GeV}/c\), it was produced in the current fragmentation region, as defined by the c.m. Feynman variable \(x_F\approx z-(E^h_T)^2/(zW^2)>0\), and it was associated with a cluster in one of the hadron calorimeters HCAL1 or HCAL2 with an energy deposit greater than 5 GeV in HCAL1 or 7 GeV in HCAL2. The efficiencies of the calorimeters above these energies are close to 100%. The energy of hadrons extended up to 120 GeV in the former and up to 140 GeV in the latter. All hadrons of the SIDIS candidates were included in the analysis of asymmetries. For the final selection of the SIDIS events and hadrons, the SIDIS candidates have to pass stability tests, as described in Sects. 3.3 and  3.4. The total number of hadrons in 2006 after afore-mentioned selections is \(15.6\times 10^6\) including \(8.6\times 10^6\) \(h^+\) and \(7.0\times 10^6\) \(h^-\).

To summarise, the SIDIS events and hadrons have been selected from preselected candidates requiring: 140 GeV/c \(< p_\mu<\) 180 GeV/c, \(Q^2 > 1\,(\hbox {GeV/c})^2\), 0.1 \(<y<\) 0.9, 0.01 \(<z<\) 1, \(p^h_T>\) 0.05, \(x_F > 0\), \(E_{HCAL1}>\) 5 GeV, \(E_{HCAL2}>\) 7 GeV. In order to test the influence of the stronger \(p^h_T\) and z kinematic cuts on the “integrated” azimuthal asymmetries, they were calculated summing up all selected hadrons (see Sect. 3.6) at variance with Ref. [6]. Azimuthal asymmetries as functions of the kinematic variables x, z or \(p^h_T\) were calculated in a restricted region following Refs. [6, 7], i.e. summing up hadrons within the intervals given in Table 1. The number of hadrons within these intervals is reduced by a factor of about two compared to the total number.

The distributions of selected SIDIS events as a function of \(Q^2\) and y and of charged hadrons as a function of z and \(p^h_T\) for different data samples are presented in Fig. 3.

3.3 Tests of data stability

Taking advantage of the three-cell polarised target, stability tests for the 2006 data were performed by investigating variations of events from run to run for certain observables via ratios \(R_i\), where i is the run number, using the combined information from cells U and D denoted by \((U+D)\), and that of cell M. One expects the ratio \(R_i =(U_i+D_i)/M_i\) to be independent of luminosity, close to unity and stable from run to run. In order to confirm this expectation, the ratios \(R_{i}\) per run were obtained for the following observables that are relevant to the selection of SIDIS events and hadrons: number of SIDIS events, number of tracks per SIDIS event, number of clusters in HCAL1 (HCAL2) with \(E > 5\ (7)\) GeV, average energy of clusters in HCAL1 (HCAL2), average energy of the associated clusters per event in HCAL1 (HCAL2) and average angle \(\langle \phi \rangle \). The \(R_{i}\) values as a function of the run number were fitted by constants \(\overline{R}\) for all runs. It was found that most of these \(R_{i}\) were stable within the \(\pm 3\sigma \) limits around the average values \(\overline{R}\approx 1.05\), except for some runs and one of the periods.

The stability of the measurement of the hadron azimuthal angle \(\phi \) in the range \(\pm 180^\circ \) is essential for determination of asymmetries. Distributions of \(\phi \)-values in this range were obtained for each run of data taking and average values \(\langle \phi \rangle _i\) per run determined. The distribution of \(\langle \phi \rangle _i\) had a Gaussian shape around the mean value equal to zero for almost all runs.

3.4 Acceptance-cancelling method for calculations of cross section ratios

The modified “acceptance-cancelling” double ratio method was used to calculate the ratios of the SIDIS cross sections for positive and negative target polarisations denoted as \(\sigma _+/\sigma _-\) and the 2006 asymmetries. In this Paper, the modified double ratio method is applied in three forms: first, in a form of “acceptance-cancelling”, second, in a form of “cross section-cancelling” and, third (see Sect. 5), again in the “acceptance-cancelling” form to test the hadron yield stability.

In order to cancel acceptances, the method utilises double ratios, i.e. the product of two ratios of events. For the three-cell target the method was modified as follows. The target cell M was artificially divided in two sub-cells M1 and M2, each 28 cm long, and two pairs of cells (U and M1) and (M2 and D) are considered below. The cells in each pair have equal lengths, i.e. equal densities of deuterons, but opposite polarisations p (\(+\) or −) at a given solenoid field direction f (\(+\) or −). For each pair of cells at a given f, one can construct the double ratio using the number of selected SIDIS events or hadrons. These numbers obtained from the cell i and denoted as \(N^i_{pf}\) are usually expressed via a product of a cell luminosity (\(L^i_f\)) given by the beam intensity times the target cell material density a target cell acceptance (\(A^i_f\)) and the corresponding cross section (\(\sigma _p\)): \(N^i_{pf} = L^i_f\times A^i_f\times \sigma _p\), i.e. luminosity, acceptance and cross section are folded in the number of events. Taking this relation into account as well as the COMPASS procedure of measurements divided in two groups of runs G1 and G2, one can construct for each pair of the target cells the “acceptance-cancelling” double ratio of events (hadrons) that provides a way to unfold the \((\sigma _+/\sigma _-)^2\). Particularly, for the polarisation settings at \(f = +\), the two double ratios of event numbers constructed for the (UM1) and for the (M2, D) pairs have the forms given by Eq. (5), where events for the first (second) ratio of each pair are taken from the G1 (G2) runs.

$$\begin{aligned}&\left[ \frac{N^U_{{+}+}}{N^{M1}_{{-}+}}\right] _\mathrm{G1} \times ~\left[ \frac{N^{M1}_{{+}+}}{N^U_{{-}+}}\right] _\mathrm{G2}= \left( \frac{\sigma _+}{\sigma _-}\right) ^2_1, \qquad \nonumber \\&\left[ \frac{N^D_{{+}+}}{N^{M2}_{{-}+}}\right] _\mathrm{G1}\times ~\left[ \frac{N^{M2}_{{+}+}}{N^D_{{-}+}}\right] _\mathrm{G2}= \left( \frac{\sigma _+}{\sigma _-}\right) ^2_2. \end{aligned}$$
(5)

Substituting in the left parts of Eq. (5) the above expressions for \(N^i_{pf}\) one can see that, after “cancellations” of \(L^i_f\) and \(A^i_f\), the double ratios of events are directly related to the cross section ratios squared. Because the luminosities of cells are equal, they contribute equally to the numerators and denominators and their cancellations are expected in each of the Eq. (5) ratios. If the acceptances \(A^i_f\) of the cells are similar at the same f in the G1 and G2 groups of runs (it is subject for tests below), they are also folded equally in the corresponding number of events and “cancel” in the double ratios, i.e. it is not necessary to calculate them (see Sect. 5). In the numerator of each ratio, the number of events (hadrons) are taken from the runs with the positive target cell polarisation, while in the denominator they are taken from the runs with negative target polarisation. Thus under above conditions each ratio of events (hadrons) in the left parts the Eq. (5) is equal to the ratio \(\sigma _+/\sigma _-\), which is known to be close to unity (it is subject for test below). Hence, each double ratio in Eq. (5), which is equal to \((\sigma _+/\sigma _-)^2\), also have to be equal within statistical uncertainties and is expected to be close to unity. The stability of acceptances during the G1 and G2 runs have been checked using the cross sections cancelling double ratios of events (hadrons) in the forms similar to ones given in Eqs. (11, 12) of Ref. [6] which are related to the ratios of acceptances.

Similarly, at \(f=-\) the two double ratios constructed for the same target pairs are:

$$\begin{aligned}&\left[ \frac{N^U_{{+}-}}{N^{M1}_{{-}-}}\right] _\mathrm{G1} \times ~\left[ \frac{N^{M1}_{{+}-}}{N^U_{{-}-}}\right] _\mathrm{G2}= \left( \frac{\sigma _+}{\sigma _-}\right) ^2_3,\qquad \nonumber \\&\left[ \frac{N^D_{{+}-}}{N^{M2}_{{-}-}}\right] _\mathrm{G1} \times ~\left[ \frac{N^{M2}_{{+}-}}{N^D_{{-}-}}\right] _\mathrm{G2}= \left( \frac{\sigma _+}{\sigma _-}\right) ^2_4. \end{aligned}$$
(6)

Because no requirements except time stabilities are imposed on the data, the Eqs. (5, 6) are valid and can be used for calculations of either the cross section ratios in the restricted kinematic regions or in the whole available region (see Sect. 3.6).

Thus each of four double ratios of the events (hadrons) in Eqs. (5, 6) calculated with SIDIS events are related to the squared ratio of the SIDIS cross sections for positive and negative target polarisations determined with a part of the data. When statistically averaged, they can be used to calculate asymmetries with the whole data provided that (i) acceptances in the G1 and G2 periods are indeed stable and equal, (ii) the values of the double ratios calculated for SIDIS events and for hadrons with polarisation settings at \(f=+\) and \(f=-\) are stable and equal within statistical uncertainties. These requirements were checked with SIDIS event candidates and hadrons and final selections of them were determined.

Altogether, the stability tests have shown that (i) acceptances are stable and equal during the G1 and G2 groups of runs, (ii) the double ratios in Eqs. (5, 6) calculated with SIDIS events or hadrons are stable over the data taking periods and contained inside the \(\pm 3\sigma \) corridors around the average values which are close to unity. In order to be accepted for analysis the average value per a given run of the acceptances, the angles \(\langle \phi \rangle _i\) and the double ratio values defined by Eqs. (5, 6) have to be within \(\pm 3\sigma \) limits of the corresponding mean value for all runs. Otherwise the run was rejected. The rejected runs contained about 10% of 2006 hadrons.

Fig. 4
figure 4

The 2006 “integrated” asymmetries a as functions of the azimuthal angle \(\phi \). Curves show the corresponding fits

Table 2 The \(h^+\) and \(h^-\) modulation amplitudes of the integrated azimuthal asymmetries obtained from the statistically combined amplitudes of the 2002–2006 (left), those of 2002–2004 (middle) and those of 2006 (right)
Fig. 5
figure 5

The values of modulation amplitudes a together with their uncertainties obtained from the fits of the integrated asymmetries \(a_{h^\pm }(\phi )\) by the function from Eq. (4) separately for the data of 2002, 2003, 2004 and 2006 as well as statistically combined modulation amplitudes for four years denoted by AV (see Sect. 4.1)

3.5 Extraction of azimuthal asymmetries in hadron production

For the extraction of the azimuthal asymmetries \(a_{h^\pm }(\phi )\) off the cross section ratios, the distributions of the charged hadrons \(h^+\) and \(h^-\) were separately analysed as a function of the azimuthal angle \(\phi \) in the region from \(-180^\circ \) to \(+180^\circ \) divided into 10 \(\phi \)-bins. For both \(h^+\) and \(h^-\), the double ratios of the hadrons defined by Eqs. (5, 6), \((\sigma _+/\sigma _-)^2_k(\phi )\)\(k=1,...4\), were calculated and combined as follows:

$$\begin{aligned}&\left( \frac{\sigma _+}{\sigma _-}\right) _{h^\pm }^2(\phi )\nonumber \\&\quad = \left[ \left( \frac{\sigma _+}{\sigma _-}\right) ^2_1\oplus \left( \frac{\sigma _+}{\sigma _-}\right) ^2_2 \oplus \left( \frac{\sigma _+}{\sigma _-}\right) ^2_3 \oplus \left( \frac{\sigma _+}{\sigma _-}\right) ^2_4\right] _{h^\pm } (\phi ) \nonumber \\&\quad \cong 1+a_{h^\pm }(\phi ) \sum _{k}\left[ \sum _{i,f,p\in k}{\mathscr {P}}^i_{pf}(x)\right] _{h^\pm }\cdot W_k~, \end{aligned}$$
(7)

where the symbol \(\oplus \) means statistically weighted averaging. As it was shown in Refs. [6, 7], in first approximation, the squared ratios of cross sections \((\sigma _+/\sigma _-)^2_{h^\pm }(\phi )\) are related to the asymmetries \(a_{h^\pm }(\phi )\) multiplied by polarisation terms. For each hadron charge, the polarisation term is given by the sum of the \({\mathscr {P}}^i_{pf}(x)\) values, each of them being the product of target cell polarisations \(|P^i_{pf}|\) and dilution factor \(\hbox {f}^{\,i}(x)\), as defined in Refs. [6,7,8], where i, p and f are those used to calculate the ratio \((\sigma _+/\sigma _-)^2_k(\phi )\), i.e. four polarisation values at each k. The weight \(W_k\) is equal to the ratio of the number of hadrons, \(N_k\), to the total number of hadrons, \(N_\mathrm{tot}\). Therefore, the \(a_{h^\pm }(\phi )\), referred to as single-hadron asymmetries, are:

$$\begin{aligned} a_{h^\pm }(\phi )\cong {\large \frac{(\frac{\sigma _+}{\sigma _-})_{h^\pm }^2(\phi )-1}{\sum \nolimits _{k} [\sum \nolimits _{~i,f,p\in k}{\mathscr {P}}^i_{pf}(x) ]_{h^\pm }\!\cdot W_k}}. \end{aligned}$$
(8)

3.6 The 2006 asymmetries

Following Sect. 3.2, the asymmetries \(a_{h^\pm }(\phi )\) were calculated in this Paper (i) as the “integrated” asymmetries using the total number of \(h^+\) or \(h^-\), and (ii) as the asymmetries vs. one of kinematic variables x, \(p^h_T\) or z disregarding the others and using the numbers of \(h^+\) or \(h^-\) within the intervals defined in Table 1. In each case, the asymmetries were fitted by the function from Eq. (4) using the standard least-square-method and extracting all asymmetry modulation amplitudes simultaneously.

Integrated asymmetries. These asymmetries for the 2006 data as a function of the azimuthal angle \(\phi \) are shown in Fig. 4 together with results of the fits given in Table 2.

In order to compare the 2006 integrated asymmetries to those of the 2002, 2003 and 2004, the latter ones were recalculated using the total number of hadrons. The results of the fit for the 2006 integrated asymmetries together with those of the 2002, 2003 and 2004 calculated similarly are shown in Fig. 5. The modulation amplitudes obtained for each year are in agreement with one another, as confirmed by the compatibility tests (see Sect. 5).

In order to check the impact of the disregarded modulations, which could have appeared from SIDIS partial cross sections \(\hbox {d}\sigma _{00}\) and \(\hbox {d}\sigma _{L0}\), on the modulation amplitudes in Eq. (4), we have performed fits of the 2006 integrated asymmetries by a new fitting function which contains a numerator and denominator. In the numerator we have used the same modulations as in Eq. (4), but in the denominator we included the disregarded modulation with average amplitudes determined in Ref. [10]. For the asymmetry \(a_{h^-}(\phi )\), the fitting function is as follows:

$$\begin{aligned} a_{h^-}(\phi )=\frac{a_{h^-}^{0} + a_{h^- }^{\sin \phi }\sin \phi + a_{h^-}^{\sin 2\phi }\sin 2\phi + a_{h^-}^{\sin 3\phi }\sin 3\phi + a_{h^-}^{\cos \phi }\cos \phi }{1+0.0412\cos 2\phi + 0.0552\cos \phi +0.0008\sin \phi }.\nonumber \\ \end{aligned}$$
(9)

By comparing results of this fit with results of the standard fit of the 2006 data shown in Fig. 5, it was found that the differences between values of modulation amplitudes in the numerator are smaller than 1% of the fit uncertainties. Similar results are obtained for \(a_{h^+}(\phi )\) replacing amplitudes in the denominator of Eq. (9) by corresponding values from Ref. [10]. Thus the contributions of the disregarded modulations to the integrated asymmetries in Eq. (4) are indeed negligible.

Asymmetries as functions of kinematic variables. The 2006 modulation amplitudes as functions of kinematic variables were compared to those from the combined 2002 – 2004 data and found to be in agreement within uncertainties of the fits. They are used for calculations of the combined 2002– 2006 modulation amplitudes (see Sect. 4.2).

4 Azimuthal asymmetries for the combined deuteron data

4.1 Integrated asymmetries

For the integrated asymmetries, the values of the combined 2002–2006 modulation amplitudes, which are shown in Fig. 5 and denoted by AV, were obtained using a statistical combination of four corresponding amplitudes. They are given in Table 2. The combined 2002–2004 modulation amplitudes, calculated as for 2006, are also shown in Table 2 in order to allow comparison to those of Refs. [6, 7], which were calculated in the restricted kinematic region. As expected for the iso-scalar deuteron target, consistent results are obtained for the \(\phi \)-independent terms \(a_{h^+}^{0}\) and \(a_{h^-}^{0}\). All \(\phi \)-modulation amplitudes are consistent with zero within uncertainties of fits. Comparing results for the 2002–2006 combined data presented in Table 2 to the results for 2002–2004 and to the results of Ref. [6], one can see that they are in agreement between themselves within the quoted uncertainties. This indicates that the integrated asymmetries with and without kinematic cuts of Table 1 are consistent, i.e. these cuts reduce statistics but do not change the values of the asymmetries within experimental uncertainties. Due to increased statistics of each year, the statistical uncertainties of the combined 2002–2006 amplitudes are reduced by a factor of about 1/1.6 compared to those of Ref. [6].

4.2 Asymmetries as functions of kinematic variables

The final 2002–2006 results on the modulation amplitudes of asymmetries \(a_{h^\pm }(\phi )\) calculated as the function of one of the variables \(x,\ z\) and \(p^{h}_{T}\) while disregarding the others were obtained from the statistically averaged 2002, 2003, 2004 and 2006 modulation amplitudes. The results are presented in Fig. 6.

Except for the \(a^{0}_{h^\pm }(x)\), all amplitudes when fitted by constants are found to be consistent with zero within statistical uncertainties (\(\chi ^2/NDF\simeq 1\)). As expected, the \(a_{h^\pm }^0(x)\) for deuteron target have the same x-dependence for positive and negative hadrons. Additionally, the x-dependence of the \(a^0_{h^\pm }(x)/D_0(x,y)\) values are presented in Fig. 7, where \(D_0(x,y)\) is the virtual-photon depolarisation factor for each x interval multiplied by the average beam polarisation \(|P_\mu |\), as defined in Refs. [6, 7]. If the \(a_{h^\pm }^0(x)\) represent the main contributions to the asymmetries of Eq. (3), the values of \(a_{h^\pm }^0(x)/D_{0}(x,y)\) by definition (see e.g. Ref. [25]) should be equal to the asymmetries \(A^{h^\pm }_{1d}(x)\). Within experimental uncertainties, there is a good agreement between our data on \(a_{h^\pm }^0(x)/D_{0}(x,y)\) and the data of Ref. [26] on \(A^{h^\pm }_{1d}(x)\), which confirms the correctness of the results on the asymmetries calculated by the modified acceptance-cancelling method. The values of \(A^{h^\pm }_{1d}(x)\) were obtained with the 2002–2004 data. A similar x-dependence was also observed with 2002–2006 data for the asymmetries \(A^{\pi ^\pm }_{1d}(x)\) and \(A^{K^+}_{1d}(x)\) obtained with the identified hadrons (see Ref. [21]).

Fig. 6
figure 6

The modulation amplitudes a of the \(h^+\) and \(h^-\) azimuthal asymmetries as the function of \(x\ ,z\) and \(p^h_T\) obtained from the combined 2002–2006 data on the muon SIDIS off longitudinally polarised deuterons. Only uncertainties of fits are shown

5 Systematic uncertainties

The global compatibility test of the results on the asymmetries \(a_{h\pm }(\phi )\), that were obtained separately for 2002, 2003, 2004 and 2006 years, was performed by building the pull distribution: \({pull}_i=(a_i-\langle a\rangle )\times |\sigma _{a_i}^2- \sigma _{\langle a\rangle }^2|^{-1/2}\), where \(a_{i}\) is the asymmetry for a given year, hadron charge and kinematic interval, \(\langle a\rangle \) is the corresponding weighted mean value over four years and \(\sigma \) denotes the corresponding standard deviation. The pull distribution had in total 750 entries compared to 540 for 2002–2004 years. The asymmetries reported in this Paper, in principle, could have the non-cancelled-acceptance or luminosity time-dependent effects folded in the event numbers of one of the Eqs. (5, 6)-ratio and, consequently, in one or several values of \(a_i\) distorting the pull distribution. In the absence of such effects, as expected, the pull distribution follows the Gaussian distribution with the mean value consistent with zero (\(-\,0.038\pm 0. 033\) in our case) and \(\sigma \) with unity (\(0.978\pm 0.024\)). This indicates that no significant time-dependent systematic effects are present in the data on asymmetries, i.e. systematic uncertainties in values of \(a_i\) are smaller than statistical ones.

Quantitative measures of possible systematic effects have been obtained by estimating additive and multiplicative uncertainties. Main contributions to possible additive systematic uncertainties could come from the instabilities of the hadron yields. The \(\phi \)-stability of hadron yields in the 2006 data was checked following the procedure described in Refs. [6, 7]. For this purpose, the double ratios of hadron numbers as a function of the azimuthal angle \(\phi \) for different polarisation settings at the field f during the G1 and G2 runs were calculated as follows:

$$\begin{aligned} f= & {} +:~~ F_+(\phi )=\frac{N_{{+}+}^U + N_{{+}+}^{D}}{N_{{-}+}^{M}}\cdot \frac{N_{{+}+}^M}{N_{{-}+}^U+N_{{-}+}^D}~, \nonumber \\ f= & {} -:~~ F_-(\phi )=\frac{N_{{+}-}^U + N_{{+}-}^{D}}{N_{{-}-}^{M}} \cdot \frac{N_{{+}-}^M}{N_{{-}-}^U+N_{{-}-}^D}~. \end{aligned}$$
(10)
Fig. 7
figure 7

The x-dependence of the values \(a_{h^\pm }^0(x)/D_{0}(x,y)\) for 2002–2006 data compared to the data of Ref. [26] on the asymmetries \(A^{h^\pm }_{1d}(x)\)

Fig. 8
figure 8

The \(\phi \)-dependence of the weighted sums of double ratios \(F(\phi )\) for 2006 data: \(h^-\), \(h^+\) and \(h^++ h^-\). The solid (red) lines represent the results of fits by constants

Here,  \(N_{{p}f}^{i}\) is the number of hadrons per \(\phi \)-bin from target cell i with polarisation p and field f, as explained in Sect. 3.4. The ratios given by Eqs. (10) are modifications of ratios used in Refs. [6, 7] for the case of two target cells. In the above double ratios, we expect the acceptance and luminosity cancellations and, as a result, the \(\phi \)-stability of the hadron yields. If unstable, they could indicate possible systematics in the acceptance as well. The fits by constants (see Fig. 8) of the weighted sums \(F(\phi )=F_+(\phi )\oplus F_-(\phi )\) in the \(\phi \)-region from \(-180^\circ \) to \(+180^\circ \) for \(h^+\), \(h^-\) and \(h^+ +h^-\) of the 2006 data gave results consistent with unity within statistical uncertainties of the order of 0.001. This means that no \(\phi \)-instabilities and acceptance-changing (not cancelled) effects have been observed, i.e. there are no large additive systematic uncertainties in the 2006 data. The value \(\Delta a_{h^\pm }(\phi )=\pm 0.001\) was chosen as a quantitative measure of possible additive systematic uncertainties in the 2006 asymmetry measurements. It is equal to \(\pm \sigma \) of the \(h^{+}+h^{-}\) data stability test for \(F(\phi )\). The same value was obtained for the 2002–2004 data [6, 7] and hence adopted also for the combined 2002–2006 deuteron data.

Possible sources of multiplicative systematic uncertainties of the asymmetry evaluation are uncertainties in the determination of the beam and target polarisations and estimations of the dilution factor. The multiplicative systematic uncertainties of the extracted asymmetries due to uncertainties in the determination of the beam and target polarisations were estimated to be less than 5% each and those due to uncertainties of the dilution factor to be less than about 2%. When combined in quadrature, overall possible multiplicative systematic uncertainty of less than 6% was obtained.

6 Conclusions

The searches for possible azimuthal modulations in the single-hadron azimuthal asymmetries \(a_{h^\pm }(\phi )\), as a manifestation of TMD PDFs describing the nucleons in the longitudinally polarised deuteron, have been performed using all COMPASS deuteron data and the acceptance-cancelling method of analysis. For each hadron charge, beside the \(\phi \)-independent term, four possible modulations predicted by theory (\(\sin \phi ,\ \sin 2\phi ,\ \sin 3\phi \) and \(\cos \phi \)) and their dependence on kinematical variables are considered. The asymmetries have been calculated both for all selected hadrons (“integrated” asymmetries) and for hadrons as functions of kinematic variables within the restricted region.

For the “integrated” asymmetries, it was found that results in the restricted range of kinematic variables are consistent with those of the wider range. In other words, the restricting of the kinematic region reduces the statistics but does not change the values of the asymmetries beyond the sensitivity of this experiment. The same result was obtained in Ref. [6] for the asymmetries as a function of z.

The \(\phi \)-independent terms \(a_{h^\pm }^0(x)\) of the asymmetries \(a_{h^\pm }(\phi )\), which are expected to originate mostly from the known helicity PDFs \(g_{1L}(x)\equiv g_1(x)\), are connected to the virtual photon asymmetry \(A^{h^\pm }_{1d}(x)=a_{h^\pm }^{0}(x)/D_0(x,y)\). There is good agreement between the COMPASS data on \(a^{0}_{h^\pm }(x)/D_0(x,y)\) and \(A^{h^\pm }_{1d}(x)\) from Refs. [21, 26] which confirms this expectation.

No statistically significant dependences of \(\phi \)-modulation amplitudes were observed as functions of x, z or \(p^h_T\) when fitted by constants. Still, there are some hints (statistically not confirmed) for a possible x-dependence of the \(\sin 2\phi \), \(\sin 3\phi \) and \(\cos \phi \) modulation amplitudes. The \(\sin 2\phi \) amplitude for \(h^{-}\) is mostly positive and rises with increasing x, while for \(h^{+}\) it is mostly negative and decreases with x. This behaviour agrees with that discussed in Refs. [19, 22, 27], if one takes into account different definitions of asymmetries by the HERMES and COMPASS Collaborations. The increase with x of the modulus of the \(\cos \phi \) amplitudes, related to the Cahn-effect [28, 29] and predicted in Ref. [30], was already visible from the 2002–2004 data [6] and persists for the combined 2002–2006 data. Hints for a possible x-dependence of \(\sin 3\phi \) modulation amplitudes are discussed in Ref. [19]. Quantitative estimates of a possible contribution of the \(\cos \phi \) modulation to the deuteron asymmetries, related to TMD PDFs \(g^\perp _L\) and \(e_L\), have been obtained in Ref. [31]. They are in agreement with our data.

Altogether, one can conclude that contributions of TMD PDFs convoluted with FFs to the azimuthal asymmetries in the cross sections of hadron production in muon SIDIS off longitudinally polarised deuterons are small. This is either due to possible cancellations of the contributions to the asymmetries by the deuteron’s up and down quarks, or/and due to the smallness of the transverse component of the target polarisation and of the suppression factor that behaves as M / Q. Some of these conclusions can be checked studying these asymmetries in muon SIDIS off longitudinally polarised protons.